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We discuss in detail a method to study transverse momentum dependent parton distribution 
functions (TMDs) using lattice QCD. To develop the formalism and to obtain first numerical results, 
we directly implement a bi-local quark-quark operator connected by a straight Wilson line, allowing 
us to study T-even, "process-independent" TMDs. Beyond results for ^-integrated TMDs and quark 
densities, we present a study of correlations in x and k±. Our calculations are based on domain 
wall valence quark propagators by the LHP collaboration calculated on top of gauge configurations 
provided by MILC with Nf — 2+1 asqtad-improved staggered sea quarks. 



I. INTRODUCTION 

The modern approach to the intrinsic quark and gluon 
structure of hadrons, in particular the nucleon, rests on 
two pillars, the generalized parton distributions (GPDs) 
[TH5] and the transverse momentum dependent distribu- 
tion functions (TMDs) 1 [SHIO]. The theoretical status 
of GPDs is fairly clear: They can be analyzed within 
the framework of collinear factorization, and have exact 
and unambiguous definitions based on off-forward hadron 
matrix elements of gauge- invariant quark and gluon oper- 
ators that are bi-local along the light cone. Transformed 
to coordinate (impact parameter, b±-) space, GPDs have 
standard interpretations as partonic probability densities 
in the longitudinal momentum fraction x and b± 
Moreover, they fully incorporate the well-known hadronic 
form factors and the PDFs, which can be obtained from 
the GPDs from integrations over x and in the forward 
limit (i.e., by integration over b±), respectively. Impor- 
tantly, at leading-twist level, all-order QCD-factorization 
theorems have been established that directly relate the 
GPDs to particular hard exclusive scattering processes 
like deeply virtual Compton scattering (DVCS) [15]. In 
this sense, the GPDs are process-independent, universal 
quantities. Moments of GPDs have been studied in lat- 
tice QCD since 2002, and for a review we refer to [13] ■ A 
calculation of GPDs performed in the same lattice frame- 
work as employed in this work has been presented re- 
cently by the LHP collaboration in Ref. |14) . 

Complementary information on the structure of 
hadrons is encoded in the TMDs. Naively, they can 
be thought of as having a probabilistic interpretation 
and describing the distribution of, e.g., quarks in a 
nucleon with respect to x and the intrinsic transverse 
momentum k± carried by the quarks, as illustrated in 
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Fig. [T] A great deal of the motivation to study TMDs 
hinges on their expected direct relation to the well- 
known "integrated" PDFs by an integration over k±. 
TMDs play a central role in h 
the description and under- 
standing of semi-inclusive 
deep inelastic scattering 
(SIDIS) processes and rela- 
ted single-spin asymmetries. 
Apart from this common 
folklore, however, one finds 
that the theoretical situa- 
tion concerning TMDs is, 
in contrast to the GPDs, 
much more challenging. In 
contrast to the GPDs, the 
framework the TMDs are 
embedded in goes beyond 

collinear factorization, and the theoretical concepts 
needed have not yet been fully developed 2 . To explain 
some of the challenges in more detail, we begin with the 
definition of a basic, momentum dependent correlation 
function, 
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FIG. 1: Illustration 
of the transverse mo- 
mentum distribution of 
quarks in the proton. 
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x-(P,S\ q(l)TU[Ci] q(0) \P,S) , (1) 
¥p(l,P,S;C) 

where \P, S) is a nucleon state of momentum P and spin 
S and r represents some Dirac matrix to be specified be- 
low 3 . The Wilson line U[Ci] is essential in order to ensure 
the gauge invariance of the expression. As usual, it can be 
represented by a path ordered exponential, see Eq. ( Al I. 
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1 also denoted as "unintegrated" PDFs. For an overview and more 
references, see also [5]. 



2 For a recent attempt in this direction see |15| 

3 For better readability, we will frequently omit the arguments q, 
P, S and C in the following. 
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FIG. 2: Illustration of the leading contribution to SIDIS. 



In a frame where the nucleon has a large momentum in +- 
direction (cf. appendix |A|) , k~ is suppressed by a factor 
~ 1/P + , and it is sufficient to consider the k~ -integrated 
correlator 

^(x,k ± ;P,S;C) = ( dk-^(k,P,S;C)\ k+=xP+ . 



(2) 

Based on its symmetry transformation properties (cf. ap- 
pendix [C]), this correlator can be parametrized in terms 
of real-valued TMDs fi, q (x, k\;C), gi, q (x, fej_; C), etc. 
[21 [lOl [16] . Concrete examples will be given in section [TT[ 
As we will see in the following, the correlator in Eq. (nT, 
and in turn the TMDs parametrizing it, will in general 
depend non-trivially on the form of the path C along 
which the quark fields at the origin and at I are con- 
nected. The question that comes to mind is if the form 
of the path is in fact uniquely determined in some way, 
or in the other extreme completely arbitrary. From a 
theoretical perspective, as long as the operator can be 
regularized and renormalized (including possible neces- 
sary modifications of the basic definition Eq. |l])), we are 
in principle allowed to consider any path we like to probe 
the internal structure of the nucleon in such a framework. 
Of strong immediate interest are of course the types of 
correlators and paths that can be directly related to ex- 
perimental observables. 

A prominent example is the SIDIS process illustrated 
inFig.[2j e.g. n(P)+j*(q) -> h(P h )+X, in a kinematical 
region where the photon virtuality is large, Q 2 = —q 2 >• 
m 2 N , and the measured transverse momentum of the pro- 
duced hadron is Ph± ~ O(Aqcd)- In this context, it 
is well known that the Wilson line U[Ci] generically rep- 
resents gluon mediated interactions of the struck quark 
with the nucleon remnants. More precisely, in pertur- 
bation theory, these final state interactions correspond 
to diagrams where arbitrarily many gluon lines are ex- 
changed, as indicated in the upper part of Fig. [2] From 
the resummed gluon exchanges (see, e.g., [17]), one ob- 
tains at tree-level a Wilson line that has the form of a sta- 
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FIG. 3: (a) Staple-shaped gauge link as in SIDIS and DY. 
(b) Straight gauge link. 



pie of infinite extent, as depicted in Fig. 3a running along 



the light-cone to infinity and back. With straight Wilson 
lines denoted by U[y,z], the staple shaped gauge link is 

given by VI [C^ 00 ^] = U[l, oov + l]U '[gov + I, oov]£/[ooi>, 0], 
where the direction v is lightlike, i>sidis = n ■ I m ~ 
portantly, it is not possible to "gauge away" effects of 
the Wilson lines by choosing, e.g., the light cone gauge 
n ■ A = 0, since the transverse part of the gauge link, de- 
pending on the gauge fields at infinity, contributes in such 
gauges [TBll21| . Furthermore, it is essential to note that 
the form of the path depends on the type of process un- 
der consideration. In particular, it turns out that in the 
Drell-Yan (DY) process, initial state interactions lead to 
a gauge link that is again staple-like but oriented in the 
opposite direction, i>dy = — ^sidiSi i-e- one finds past- 
in contrast to future-pointing Wilson lines [2D] . These 
well known observations clearly show that even in a phe- 
nomenological context, already at tree- level in perturba- 
tion theory the form of the Wilson line connecting the 
quark fields in Eq. ([I]), and therefore the structure of 
the correlation function itself, is non-unique. On the 
level of the TMDs, the different directions v for SIDIS 
and DY translate for example into a sign change of so- 
called time reversal odd TMDs such as the Sivers func- 
tion, /^(^fci;^ 00 ™)) = -f^x^l;^- 00 ^). The im- 
portant message is that the TMDs can therefore be seen 
as non-universal objects, albeit the "breaking" of uni- 
versality is exactly calculable, at least in the considered 
cases. Another way of formulating these observations is 
to consider linear combinations (the sum and difference) 
of future- and past-pointing Wilson-line operators, lead- 
ing to "T-even" and "T-odd" correlators that are sep- 
arately process- independent. The non-universality can 
then be seen in the fact that there exist two distinct 
classes of TMDs, the T-even and T-odd TMDs, which are 
based on two types of operators with fundamentally dif- 
ferent gauge link structures [35] . We note that additional, 
even more complex gauge- link structures have been found 
in the framework of tree-level analyses of 2 — > 2 hadron 
scattering processes |23j . However, a more recent study 
[2~3] argues that a generalized TMD factorization of this 
kind (see also [T71[5S]) cannot be achieved for such pro- 
cesses. The argument is based on a model calculation 
that gives an explicit example where it is impossible to 
find standard Wilson line structures that allow factoriza- 
tion. 



3 



In summary, for SIDIS and the Drell-Yan process at 
tree-level one finds a standard factorization of hard and 
soft parts, where the latter, illustrated in the lower part 
of Fig. [2j is represented by the correlator in Eq. ([!]), 
with a Wilson-line of the form shown in Fig. [3a[ This 
picture changes completely as soon as loop-corrections 
are taken into account in the lower part of Fig. [2j Al- 
ready at leading one-loop level, one finds that the light- 
like sections of the Wilson lines lead to divergences due 
to light-cone singularities in the additional gluon prop- 
agator [55]. Hence, to obtain well defined amplitudes, 
the basic definition in Eq. with the staple-like Wilson 
line along the light-cone has to be modified. Different 
improved definitions of TMDs and strategies to remove 
the divergences have been proposed and discussed in the 
literature [51 I2"7rf3"2"] . To illustrate the theoretical status 
of these issues, we briefly discuss in the following two 
different approaches. In [28], a QCD factorization theo- 
rem for SIDIS has been established at leading one-loop 
level 4 , where the vector v has been taken slightly off the 
lightcone (i.e., timelike) to regularize the light-cone di- 
vergences. This leads to an additional dependence of the 
correlators on the energy of the incoming hadron, or the 
variable ( = (2P-v) 2 jv 2 , which is described by a known 
evolution equation in certain kinematical regions. Fur- 
thermore, in order to cancel out extra soft contributions 
from the basic correlator, the definition Eq. ([I]) has to 
be modified to include appropriate vacuum expectation 
values of products of Wilson lines. An important point 
is, however, that in this approach the light-cone limit, 
v 2 —> 0, cannot be taken exactly, and that no direct rela- 
tion to the standard PDFs, e.g. through an integration 
over k±, can be established. This leads clearly to some 
tension with respect to the increasing number of phc- 
nomenological analyses and parametrizations of SIDIS 
experiments (e.g., in Ref. [33]), which on the one hand 
should be based on a QCD factorization theorem, but on 
the other hand, so far make use of the assumption that 
the involved TMDs reduce to the PDFs after integration 
over k±. 

An alternative definition of TMD-correlators has been 
worked out in Ref. [51] . It is based on an exactly light- 
like direction v and a different regularization of the light- 
cone singularities involving certain pole-prescriptions. In 
order to remove the prescription dependence at least 
at one-loop level, sections of the gauge-link path that 
run along the transverse direction to infinity, i.e., from 
(oov + 0j_) to (oot> + 00 jj and back to (oov + have 
to be explicitly taken into account . In addition, a soft 
counter term has to be included in the modified definition 
of the correlation function in Eq. Q. A clear advantage 



of this approach is that the (dimensionally regularized) 
fc_i_-integral of the TMDs defined in this way reproduces 
the standard PDFs. However, it is not known to this 
date if the TMD-correlator defined in Ref. [3T] is part 
of any QCD-factorization theorem of a physical process, 
which would be a necessary condition for any solid phc- 
nomenological analyses. 

In summary, the current situation turns out to be quite 
challenging. Finding a definition of TMDs that allows to 
relate them to the PDFs, and that at the same time is 
part of a proper factorization theorem for, e.g., SIDIS, is 
non-trivial and still a matter of ongoing research. 

In view of the issues discussed so far, and the impor- 
tance of TMDs for our understanding of hadron struc- 
ture, we propose to start a program of systematic non- 
perturbative studies of the relevant correlation functions 
in the framework of lattice QCD, in addition to the on- 
going perturbative investigations. Keeping in mind that 
the lattice discretization of QCD represents a manifestly 
gauge-invariant scheme with build-in cutoff, and that the 
non-perturbative evaluation of the path integrals doesn't 
require a fixing of the gauge (which in the perturbative 
analyses contributes substantially to the difficulties), the 
lattice approach has the potential to provide new insights 
into the general properties of possible TMD-correlators 
from a completely different perspective. The long term 
plan is to perform non-perturbative studies of matrix el- 
ements of manifestly non-local operators with different 
gauge-link structures, of potentially relevant soft factors 
(vacuum expectation values of Wilson-lines and -loops) , 
and to get quantitative information from first principles 
about the x-and fcj_-dependences of the TMDs. 

The direct implementation of non-local operators like 
q(l) T U[C{\ q(0) on the lattice is still a novelty. There- 
fore, our first steps will be based on simplified operator 
structures, allowing us to establish the basic ideas, for- 
malism and methodology, and to perform first studies 
of lattice related issues like the renormalization of po- 
tential power-divergences of the Wilson-lines and certain 
discretization effects. Specifically, taking into account 
the fact that there is no straightforward way to realize 
lightlike gauge links on the lattice, we have performed 
first investigations with a simple path geometry: We em- 
ploy a direct, straight Wilson line U[Cf w ] = U[l,0], see 
Fig. 3b The straight Wilson line ("sW") is a process- 



4 The validity with respect to higher order corrections is still under 
debate. 

5 In contrast to the covariant gauge used in Ref. |28| , the trans- 
verse sections at cod cannot be neglected in the light-cone gauge 
that was employed in this case. 



independent choice that serves us here as a starting point 
for exploratory calculations. Note that time reversal odd 
TMDs vanish by symmetry for straight Wilson lines, e.g., 
f^ T (x,k ± ;C sW ) = 0. Although our TMDs defined in 
this way are thus not directly related to those defined 
and used in the literature and for the description of, e.g., 
SIDIS, they still can be seen as being elements of the 
general class of "process-independent, T-even" TMDs, as 
discussed above. Although being preliminary, our com- 
putations therefore provide some semi-quantitative in- 
formation about this class of TMDs, in particular with 
respect to their signs and (relative) sizes. First numeri- 
cal results have already been presented by us in Ref. |34) , 
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where we observed clear signals for several TMDs, corre- 
sponding to sizable correlations in k± and the quark and 
nucleon spins, s and S 1 , leading to visibly deformed den- 
sities of (polarized) up- and down-quarks in a (polarized) 
nucleon. Here, we give a more detailed description of our 
techniques, and discuss critical issues as well as possible 
improvements and extensions. 



II. PARAMETRIZATION IN TERMS OF TMDS 
AND INVARIANT AMPLITUDES 

We now come back to the parametrization of the 
k~ -integrated correlator in Eq. ([2| in terms of TMDs. 
Following the common conventions in the literature 
[3 13 EH US]; we decompose the correlator for T = 
7+,7+7 5 , ia l+ j 5 into the leading twist-2 TMDs as fol- 
lows: 



^^\x,k ± ) = ^{^h± 

P + TUN 



<f>^(x,k ± ) = h- 
^ + ^\x,k x ) = A gi + 



£jj kj Sj 
m N 

k±S x 



m N 



9lT 



(x, k ± ) =S t h 1+ hyp 



m N 



2m 2 N 
m N 



(3) 
(4) 

(5) 



Here i, j = 1,2 are indices denoting transverse directions. 
The TMDs in square brackets are odd under time rever- 
sal and absent for our choice of a straight Wilson line. 
For other Dirac structures T, the correlator $I r l (a;, fcj_) 
is suppressed by factors / P + or (m» / 'P + ) 2 , corre- 
sponding to contributions of higher twist-3 and twist-4, 
respectively. The parametrizations of the twist-3 corre- 
lators are given by [SJ [TU1 [TB] 
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Cijh 



*^%,k x ) = ^ \Ah L + k ^ S ±h T 
P+ m N 



(10) 
(11) 



where square brackets around pairs of indices denote an- 
tisymmetrization, a^b^ = a^V — a u b^. Naively, one 
might ask how the TMDs defined in Eqns. to fllT) ), 
that are classified according to twist and part of an ex- 
pansion of correlators in mjv/P + with large P + , can 
ever be accessed in lattice QCD simulations, where the 
nucleon is at rest or has only a small non-zero three- 
momentum. A first step towards the resolution of this po- 
tential contradiction is a frame independent parametriza- 
tion of $ [ P(1,P,S;C) on the right hand side in Eq. Q 
in terms of Lorentz- invariant amplitudes Ai(l 2 ,l-P). As 
will be explained in the following sections, the non-local 
operator technique allows us to evaluate the Z-dependent 
matrix element $ q (I, P, S;C) directly on the lattice. 

Analogous to the procedure outlined in Ref. [55], we 
write down all Lorentz-covariant structures compatible 
with the properties of <£>q (l,P, S;C) under symmetry 
transformations, see appendix[C] For straight gauge links 
C sW , we obtain: 



$W =2m N A 1 , 
$^ 5 i = , 

$[t"] = 2 P^ A 2 + 2i m N 2 V 1 A 3 , 
$rW] = -2m N S"A 6 -2im N Pi J -{l-S)A 7 
+ 2m N 3 F(l- S)A 8 , 
^•"V] _ 2 pfrgv] 2 g + 2i m 2 N l^S u] A w 

+ 2m? N &P v \l- S)A X1 . (12) 

The structures above can be obtained by replacing k 
by im 2 N l in the corresponding structures for the time- 
reversal-even amplitudes A l in Ref. QI]. 6 The rep- 
resentation in terms of the Ai(l 2 ,l-P) is a more con- 
venient choice for our purposes than the conventional 
parametrization using momentum dependent amplitudes 
Ai(k 2 ,k-P). The (I 2 , Z-P)-dependent representation will 
also be advantageous for the discussion of correlation s in 
the x- and fc^-dependence of the TMDs, see section 
The amplitudes A t are complex-valued and fulfill 



VI 



Ai(l 2 ,l-P) = A z (l 2 ,-l-P) 



(13) 



This property follows from hermiticity and is analogous 
to the constraint that the TMDs and the conventional 



(9) 



6 We adjust our sign conventions for Ag, Aiq and An as well as 
the linear combination Ag m = Ag — ^m 2 N l 2 A\\ with respect to 
previous work 134 , 36 in favor of this simple correspondence. 
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amplitudes Ai(k 2 , k-P) are real. Notice that there is in 
general no one-to-one correspondence between an indi- 
vidual Ai(l 2 ,l-P) and the Fourier-transform of the anal- 
ogous Ai(k 2 ,k-P). For example, A 8 contributes to A 6 , 
A-j and As (following the conventions of Ref. [TS]). 

Clearly, the momentum dependent amplitudes 
Ai{k 2 , k-P), as well as our invariant complex amplitudes 
Ai(l 2 ,l-P), contain information about all leading and 
higher twist contributions (for the given choice of the 
Wilson-line path). To see how the TMDs of different 
twist can be obtained from the invariant amplitudes, 
we first note that combining the definitions ([!]) and ([2]), 
the k~ -integral in Eq. (p]) translates into the constraint 
1+ = 0. Using r = {IP)/P + for 1+ = 0, we obtain 



d 2 l 



x ±^(l,P,S) 



(2tt) 2 



-L il±k 1 



(14) 



Inserting the structures in Eq. ( [12] ), the angular part of 
the Zj_-integral can be performed. Due to the restriction 
to l + = 0, the remaining radial integral can be rewrit- 
ten as an integral over I 2 = For the following dis- 
cussions, it is therefore useful to abbreviate the Fourier- 
transform of amplitudes as 



A, 




-L e ~ix(l-P)+il ± -k,j 



P\k±\) Mi 2 , i-p) 



(15) 



where Jo is a Bessel function. Notice that x -f-> (l-P) 
and k\ -f-> I 2 form pairs of conjugate variables with 
respect to the Fourier transform. Notice also that 
I 2 < in the Fourier-integral above. It turns out that 
only spacelike and lightlike quark separations I occur 
in the matrix elements needed for TMDs. In the fol- 
lowing, we shall use the abbrevation |Z| = V— I 2 - Fi- 
nally, the TMDs can be identified and extracted from 
comparisons of the parametrizations in Eqns. ([3])- 
([5]) with Equations (14 1 and (12), and turn out to 
be given by certain linear combinations of (x- and 
fe_i_-derivatives of) the Fourier-transformed amplitudes. 
Specifically, we obtain the twist- 2 TMDs from the ampli- 
tudes ^2,6,7, 9m, 10, 11 (I 2 , l-P): 

fi(x,k 2 L ) = 2^A 2 , 

gi (x, k 2 ± ) =-2jfl 6 + 2d x jf A 7 , 



gir(x, fejj = 4m 2 N d k 2^ ^ A 7 , 

hi L (x, fej_) = 4m 2 N d k 2^ (^A w + d x ^ Anj , 



hi(x,k 2 x ) = -2 
hi T (x, fejj = 8m 



■j' Ag m , 



An 



(16) 



Here Ag m = Ag — hm 2 N l 2 A\\. As an example for corre- 
sponding relations at subleading twist, we note that the 
axial-vector TMDs g' T and of twist-3 can be obtained 
from 

g' T (x, k 2 ± ) = -2 j A 6 + im 2 N d kl j A s , 
g£(x,kl)=8m%(d ki y fl s . (17) 



Eqns. (16) and (17) finally show that the specific types 



of linear combinations and (derivatives) of the involved 
amplitudes indeed allow a projection of the invariant Ai 
on TMDs of definite twist. 

To forestall potential confusion, we also note that the 
number of independent amplitudes in Eq. ( 12 ) (which 
is 9) is already lower than the total number of T-even 
TMDs of twist-2 and twist-3 TMDs in Eq. (p} and ilfh 



(which is 16), respectively, leaving aside the contributions 
of twist-4. This is a direct consequence of our choice of 
a straight Wilson-line path, i.e. the fact that no addi- 
tional structures depending on a direction vector v pfc I 
can appear in the parametrization E q. (|12[ ) . Accordingly, 
by a comparison of Eqns. ( 16 ) and ( |17[ ) :or example, it 
is possible to derive certain relations between (deriva- 
tives) of TMDs of twist-2 and twist-3 that are exact for 
our process-independent choice C — C sW . Such relations 
are similar but not identical to the so-called "Lorentz- 
invariance relations" [HI [Till , which only hold if the de- 
pendence on the direction vector of the staple-like gauge 
links, i.e. v = n in Fig. 3a is neglected. 
Integrating Eq. ( 14 ) over fcj^ , we obtain 



*M(x;P,S) = J ^^e-^W\l,P,S) 
dl 



l+=l ± =0 



-il-p-* 



2(2rr) 



x (P,S\ q(l-n)TU[d- n ] g (0) \P,S) 



(18) 



A parametrization of the above correlator yields the con- 
ventional, "integrated" PDFs. Notice that the staple 
shaped links of Fig. [3a] simplify to a simple connect- 
ing straight light-like Wilson line in the matrix element 
above, because the quark fields have no transverse sep- 
aration. Due to the perturbative tail of the correla- 
tor in Eq. ([T]) at large transverse momentum, the fcj_- 
intcgrations are formally divergent [38] and require a 
regularization. PDFs are typically introduced directly 
according to Eq. ( 18 ) based on renormalized operators. 
The divergent fc^-integral thus does not appear explic- 

itly 
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III. LATTICE CALCULATIONS 

A. The discretized non-local operator 

A first important step in the lattice calculation of 
TMDs is to find a discretized representation of the con- 
tinuum operator 



x 2 A 



O r , q [C l }(z) = q(l + z)TU[C l +z} q(z) 



(19) 



that appears in the matrix element $I r l of Eq. ([IJ. 
Note that we have introduced an overall offset z, 
which does not affect the matrix element: <F r ] = 
±{P,S\O r , q [Ct](z)\P,S) is independent of z. To im- 
plement the non-local operator Or, q [Ci)(z) on the lat- 
tice, we approximate the Wilson line U[Ci + z] between 
the quark fields by a product of connected link vari- 
ables, as illustrated in Fig. [3] and explained in the 
following. With the notation U^(x) = U(x,x + ae M ), 



uU 



U(x + ae^,x), the lattice gauge link for a lat- 



tice path C\ at = (x (n \x ( - n ' 1 \x ( - n - 2 \...,x ( - 1 \x^) along 
adjacent lattice sites x^> is 

U lat [Cl at ] = U(x^ ,x^ n -^) ■ ■ - U(x^ ,x^)U(x 



(i) x (o) 



) ■ 

(20) 

The above expression converges to the Wilson line Eq. 
( Al ) in the naive continuum limit, provided the distance 
of the points x^ l > to the continuous path Ci is guaranteed 
to be of the order of the lattice spacing, see appendix [5] 
As a whole, the lattice field combination we employ to 
probe nucleon structure, 



lat 



}(z) EE q(l + z)rU Ut [Cl at + z] q(z) , (21) 



has the same form as the continuum operator in Eq. ( 19 1 , 



except for the discretized gauge link along the lattice path 
Qiat runn j n g f r om the origin, a;' ' = 0, to x^ n ' = I. 

If / is a multiple of one of the unit vectors e M , C lat is 
a straight path that lies on one of the lattice axes. If I 
is at an oblique angle, we employ a method similar to 
the Bresenham algorithm J3U] to generate a step-like lat- 
tice path close to the continuum path, as in the example 
shown in Fig. [4] 

The renormalization of the lattice operators and fur- 
ther properties of the gauge link will be discussed in sec- 
tion [TTrrJl [TVBl and |r7Cl below. 



B. Lattice correlation functions 

Using the discretized non-local operator of the previous 
section, we extract the invariant amplitudes Ai(l 2 ,l-P) 
from lattice three-point correlation functions correspond- 
ing to the matrix elements $I r l. A typical lattice three- 
point-function with a non-local operator insertion at Eu- 
clidean time t is illustrated in Fig. [5j where the nucleon 
source and sink are placed at t SIC and t sn k, respectively. 

The evaluation of three-point functions follows stan- 
dard techniques 40-42 which we review very briefly in 



1. 



FIG. 4: Example of a step-like link path: The straight gauge 
link in the continuum with I = (6, 3, 0) (dashed line) is rep- 
resented as a product of link variables (7 M in the directions 
fi= 1,2,1,1,2,1,1,2,1. 



the following. Only the operators Op ? [Cj I ] we use to 
probe the nucleon and the way we interpret the results 
are specific to our task. The purpose of the source and 
the sink is to create and annihilate states with the quan- 
tum numbers of the nucleon. The nucleon sink has the 
form 



B a {t,P) 




e iP x e abc x 



*) [ U I( X ^) r diq d c (x,t) 



(22) 



where a, b, c are color indices, a is a Dirac index, T di( i = 
747275(1 + 74) and P is the three-momentum of the nu- 
cleon. An analogous expression B a (t,P) acts as a nu- 
cleon source. To increase the overlap with the nucleon, 
the quark fields u and d that enter Eq. p2| are smeared 
as described in Ref. [41]. We introduce the two-point 
function by 



C 2 Pt(p) ^J2r 2 £{B a (t snk ,P)B p (t s 

j3a 



and the three-point function for a general operator O is 
given by 



C 3pt [0 lat ](Pj T) 1 £ £ r 3pt {Ba{ts ^ p) 
L z f3a 

O ut {z,r) B p {t STC ,P)) . 



(23) 



where (■••) = J V[q, q,U] ■ ■ ■ exp(— 5" lat ) denotes an ex- 
pectation value defined by the lattice path integral, and 
where r 3pt is a Dirac matrix projecting out the desired 
parity and spin polarization of the baryon. 

In order to ensure that the transfer matrix formalism 
enables us to rewrite our three-point function in terms 
of a matrix element (N(P, S') \ lat \N(P,S)), we limit 
ourselves to operators Op at g [C; lat ](;z, r) that do not extend 
in the Euclidean time direction, i.e., the link path is re- 
stricted to the spatial plane at r, and Z 4 = 1° = 0. As 
explained in section [TTJ our selection of vectors / and P 
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on the lattice does not need to correspond to the large 
momentum frame usually chosen to introduce TMDs in 
the context of scattering processes. Relevant for the cal- 
culation of the TMDs from the amplitudes Ai(l 2 ,l-P) 
are only the Lorentz-invariant quantities formed by the 
Minkowski four-vectors I and P, which are in the lattice 



frame given by I 2 



-I 2 



III and IP 



IP. 



Consequently, we will only be able to evaluate the am- 
plitudes Ai(l 2 ,l-P) in the range 



<0, 



(24) 



where P is the chosen nucleon momentum on the lattice. 

The transfer matrix formalism shows that the lat- 
tice correlation functions decay exponentially in the Eu- 
clidean time and the energies of the contributing states. 
If the operator position r is far enough away from source 
t src and sink i sn k, the three-point function is there- 
fore dominated by contributions proportional to nucleon 
ground state matrix elements (N(P,S')\ lat \N(P,S)). 
The proportionality factors (e.g., overlaps of nucleon 
source and sink with the nucleon state), the exponen- 
tial time dependence, as well as part of the statistical 
noise cancel in the ratio with the two-point function 



R[0 lat }{P,T) 



= C3pt[ lat](p )T ) 



(25) 



If t src and i sn k are far enough apart, we observe a 
"plateau" in a region where the ground state dominates, 
such that R[O la,t ](P , t) is independent of t: 

RiO^iP.r) |T ' Wl ' |T ' t3 " kl>>A£ ' 1 > i?[0 lat ](P) , (26) 



R[0](P) = Y^ 



u(p, s) r 3 ?* u{p, S') 



%s ', 2E P tr D {r 2 Pt {-if + m N )} 
x (N(P,S')\ O \N(P,S)) , (27) 

where AE = E' — E is the difference between the ener- 
gies of the ground state and the first excited state, and 
U(P, S) is the Dirac spinor of a nucleon. For an appro- 
priately renormalized lattice operator O^, we identify 
this plateau value with the correspondingly renormalized 
continuum expression: 



RHf n }{P) ^\ R[O rcn ]{P) 



(28) 



Thus we finally gain access to the desired continuum 
matrix elements (N(P,S')\ O rcn \N(P,S)). With Equa- 



tion (27) for R[O r r ™[Ci]]{P) and inserting (for the case of 



straight gauge paths Ci) our parametrization Eq. (12 1, 



we can parametrize the plateau values in ter ms of the 
amplitudes Ai 1 as given explicitly in Table VI in the ap- 
pendix. 

We now discuss the strategy for evaluating the three- 
point function C 3pt [0^i[Cl at }]{T,P). The average over 



all offsets z in Eq. |23j) increases statistics and allows 
us to exploit translation invariance in favor of a fixed 



nucleon 
source 
(fixed position) 



nucleon 
sink 

(fixed momentum) 




time 



FIG. 5: Schematic diagram of a nucleon three-point function 
on the lattice, here for an operator probing d-quarks. 



source location. Integrating out fermions analytically, 
pairs of quark field variables u, u and d, d combine into 
lattice quark propagators, which we depict as connecting 
lines between the quark variables in Fig. [5j Lattice quark 
propagators are numerically obtained by inversion of the 
lattice Dirac operator and describe the propagation of a 
valence quark in a gauge field background, i.e., effects 
of gluons and sea quarks are included. In principle, all 
possible contractions of pairs u, u and d, d into propaga- 
tors must be taken into account. In Figure [5j a second 
diagram, resulting from the permutation of it-quarks, is 
indicated with dashed lines. In practice, however, we 
neglect here the computationally demanding so called 
disconnected contributions, where the quark variables of 
Op at ? [C/] contract with each other internally to form a 
closed quark loop. Disconnected contributions cancel ex- 
actly in the isovector case, i.e., for Oj?^, 



— /nlat 
-d — u T,u 



u r,d- 



For the numerical calculation, we employ the sequen- 
tial source technique [33], which permits us to evaluate 
the three-point function for arbitrary gauge link paths 
using the same given set of point-to-all type lattice prop- 
agators. As indicated by the curved gray envelope in 
Fig. [5j three of the quark propagators in the diagram 
can be combined into a single "sequential propagator" , 
which can be calculated for fixed (x src , i src ), i sn k and P 
using a secondary inversion, and which can be used like 
a backward point-to-all lattice propagator. Finally, the 
three-point function is evaluated by forming a product of 
a forward propagator, the link variables and the sequen- 
tial propagator. 
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Simulation parameters and computational 
details 



For the purpose of our proof-of-concept calculations, 
we have chosen existing ensembles and propagators at 
intermediate pion masses that have already been suc- 
cessfully used in the determination of GPDs [42]. The 
gauge configurations have been generated by the MILC 
collaboration [44-46]. They feature 2+1 dynamical, im- 
proved staggered quarks, with the strange quark mass 
fixed approximately to the physical value. Employing 
the "coarse" MILC gauge configurations (a ps 0.12 fm), 
the LHP collaboration has calculated propagators using 
a domain wall fermion action, where the pion mass has 
been adjusted to the Goldstone pion mass of the under- 
lying staggered lattice [32]. The computationally more 
expensive domain wall action for the valence quarks ex- 
hibits a lattice chiral symmetry, which is in particular ad- 
vantageous with respect to the operator renormalization. 
Essential ensemble parameters, together with the pion 
mass determined using the domain wall propagators, are 
given in Table [I] The MILC collaboration has chosen the 
strange quark masses m s to correspond roughly to the 
physical value. For our scaling study in Sec. |IV C[ we 
take advantage of "fine-04", "superfine-04" and "extra- 
coarse-04" gauge configurations that have become avail- 
able from the MILC collaboration recently. The ensem- 
bles listed in the last four lines of Table |T] all have the 
same ratio TO Ui( j/m s = 0.4, placing them approximately 
on a line of constant physics, i.e., they feature similar 
pion and kaon masses. In order to determine the lattice 
spacing in a uniform way for all six ensembles in the ta- 
ble, we have taken the updated, "smoothed" values rr/a 
of Ref. [47], and r\ = 0.3133(26) fm from the recent anal- 
ysis Ref. 0H]. 7 

To reduce computational costs for the production of 
propagators further, the coarse lattice gauge configura- 
tions have been chopped into two halves of temporal ex- 
tent T/2 = 32. Only every sixth trajectory, and alternat- 
ing temporal halves have been selected, reducing auto- 
correlations to an undetectable level. Noise has been re- 
duced by application of HYP-smearing [5U] to the gauge 
configurations before the propagators have been deter- 
mined by inversion. Link smearing is an operation in 
which each link variable is replaced by a unitarized "av- 
erage" of itself and gauge links in the vicinity. In the 
case of HYP-smearing, only link variables from within 
the lattice hypercubes adjacent to the original link enter 
the average, so as to minimize the distortion of physi- 
cal properties at short distances. An important benefit 
of HYP-smearing is a reduction of the breaking of rota- 




(fm) 



FIG. 6: Coverage of the (I 2 , Z-P)-plane for our choice of link 
paths. The scale on top is in lattice units and the scale on 
the right labels the integer values accessible on the periodic 
lattice. For the conversion to physical units (scales on the left 
and bottom axes) we use L/a = 20 and a — 0.1 166 fm, i.e., 
the values listed in Table [I] for the course- 10 ensemble. Note 
that l-P is dimensionless in natural units. 



7 In contrast, Refs. |34II42| used a = 0.124 fm, as determined from 
the T spectrum on the coarse lattices 451l49|. As a result, num- 
bers in physical units, including the pion masses listed in Table 
[i] differ somewhat with respect to these previous references. 



tional symmetry, see also section [IV B| below. The prop- 
agators and sequential propagators provided by LHPC 
are of the smeared-to-point type, i.e., the quark fields 
at the source location and the nucleon sink embedded 
in the sequential propagator are smeared as described 
in Ref. [31]. Using the smeared-to-point propagator as 
input, we form a smeared-to-smeared version, in order 
to be able to compute the appropriate two-point func- 
tion with smearing both at source and sink. The sequen- 
tial propagators are available for sink momenta P = 
and P = (—1,0,0) x 2n/L. The latter corresponds to 
\P\ ps 500 McV and is the lowest non-zero momentum 
on these lattices. The source-sink separation is fixed to 
4nk - hie = 10 w 1.2 fm. 

For our analysis with lattice nucleon momentum P = 
(0,0,0), we have generated 263 different link paths Cj at . 
We remind the reader that we restrict ourselves to purely 
spatial extensions of the gauge link. The quark separa- 
tions I cover the three lattice axes up to a link length 
\l\ = 20a, three quadrants in the (ii,^)- and (£1,(3)- 
planes for \l\ < 8a and a choice of additional links with 
\l\ < 15a in the first octant. For the analysis with 
P = (—1,0,0) x 2n/L, we choose 743 further vectors 
I from the two octants with I2 > 0, £3 > such that the 
Z-P)-plane is densely covered in the range accessible 
on the lattice, see Fig. [6] 

In Figure [7] we show an example plot of the ratio 
i?[Op at 9 [C; at ]](r, P) as a function of r between t src and 
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0.05930(08)(50) 
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0.082 
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TABLE I: Lattice parameters of the MILC gauge configurations [44, 45] \$T\ used in this work. The first error quoted for a 
estimates statistical errors in r\/a, the second error originates from the uncertainty about ri in physical units. The sixth and 
seventh column list the pion and the nucleon masses as determined with the LHPC propagators with domain wall valence 
fermions [42]. The first error is statistical, the second error comes from the conversion to physical units using a as quoted in 
the table. Note that the masses quoted here in physical units differ slightly from those listed in Ref. [42], because we use 
a different scheme to fix the lattice spacing, see footnote 7. The last column lists the number of configurations used for the 
calculation of three-point functions. 



t sn k- Even for the rather long link path of Fig. [4] with 
\l\ f« 0.8 fm, the signal to noise ratio is good. We fol- 
low the strategy of Ref. 02] and take the average of the 
three data points at t - t STC — 4, 5, 6 as an estimate of 

(P) defined in Eq. 



(251. 



the plateau value R[O l ^ q [C\^ 
Potential contaminations from excited states can be ne- 
glected at our present level of accuracy. 

In order to estimate statistical uncertainties, we con- 
sistently employ the Jackknife method [5H453) . For fits 
to lattice data, we minimize for each configuration j 

2 



X 



,plP)-yr> Ay, 



(29) 



Here fi denotes the fit function evaluated at location i, 
the lattice data at this location are given by Jackknife 
samples y\ , and the Jackknife error is Aj/j. The pa- 
rameter estimates pf \ ■ ■ ■ ,Pn) thus obtained are again 
Jackknife samples. The functional form of Eq. (29 1 does 



not reflect correlations among the data points by means 
of the covariance matrix. Nevertheless, the least squares 
fit using x 2 as given above implements a consistent esti- 
mator [54 for the Jackknife samples p{ , . . . ,p„ . Hence, 
the Jackknife errors that are finally obtained for the pa- 
rameters and for functions of the parameters adequately 
include correlations. 



D. Renormalization of the non-local operators 

The renormalization properties of the continuum op- 
erator Or ig [C;] have been studied with the help of an 
auxiliary field technique ("z-field") in Refs. f55H5"8"] and 
independently in leading order perturbative QCD in Ref. 
[59] . For a smooth open path C; , the renormalized Wilson 
line has the form 



U rcD [Ci] = z- l er 5mi[Ci] u[Ci 



(30) 



where £[Ci] is the total length of the path. The length 
dependent, exponential factor corresponds to the self- 
energy of the Wilson line. The dimensionful renormaliza- 
tion constant 5m removes a divergence linear in the cutoff 
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coarse-10 lattice, 



Plateau plot for the real part of the ratio 
r) as a function of r for the HYP-smeared 
for P = and the link path C ; lat depicted 
in Fig. [Z] The plateau value R[0^l u _ d [Cl at ]](P) is extracted 
from the three encircled points and is displayed as a horizontal 
error band. 



scale (i.e., a -1 on the lattice). In dimensional regulariza- 
tion, 5m vanishes, but renormalon ambiguities appear, 
see e.g., Ref. [BO]. The renormalization factor Z~ x can 
be associated with the end points of the gauge link and 
does not appear in a Wilson loop. For a piecewise smooth 
gauge link, we would have to add an angle-dependent 
renormalization factor for each corner point. For the 
composite operator Or,g[C;], we get an additional renor- 
malization factor Z2 for the quark field renormalization 
for the quark - gauge link vertices: 



and a factor Z?^ 



u r,qM\ — ^-0 (V>«0 z 



-5me[Ci] 



(31) 



J tp,z 



Note that the renormalization constants do not de- 
pend on r. This is in contrast to the renormalization 
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of local operators of the form q(0)Tq(0), q(0)T D fj,q(Q) , 
q(0)TD^D v q(0), ... as they are used, e.g., in the cal- 
culation of moments of GPDs. The basic explanation 
for the T-independent renormalization of the non-local 
object is that the spatially separated quark fields are 
renormalized individually. However, the precise relation 
between the derivative operators q(0)r_D M Z? 1/ ■ ■ ■ q(0) and 
the non-local operator Op e ^[C;] remains to be studied fur- 
ther. The interested reader is referred to appendix [HJ 
where we rewrite the non-local lattice operator Op* [Cj] 
explicitly as a weighted sum of derivative operators. The 
main purpose of appendix [H] is to address the question 
whether and how mixing among local operators affects 
the non-local object. 

As discussed in section [IV C[ it is known how to renor- 
malize straight Wilson lines on the lattice that run along 
the lattice axes. It turns out that renormalization in this 



(a) 



case is also of the form Eq. ( 30 1 . 

Fundamental for the remainder of this work, we 
will make the assumption that the discretized operator 
Op at g [Cj at ] has the same renormalization properties as the 
continuum operator. Specifically, we will employ Eq. 
(311 to renormalize our lattice operator O^JCj 8 *]. This 



assumption relies on the physical argument that, for a 
given discretization prescription of the gauge link, the 
operator O r at g [C| at ] becomes an approximate representa- 
tion of the continuum operator Or,g[Cz] as soon as the 
length of the gauge link is large compared to the lat- 
tice spacing a. Note that the numerical values of the 
renormalization constants we obtain for given renormal- 
ization conditions depend on the lattice action used and 
on the details of implementation of the discretized op- 
erator. Numerical checks of these assumptions and the 
non-perturbative methods that are employed to deter- 
mine the renormalization constants for given lattice ac- 
tion, lattice spacing and renormalization conditions will 
be discussed in sections |IVB| and |IV C| We point out 
that more detailed work on the renormalization of the 
general, step-like non-local lattice operator could bene- 
fit from the method of constructing symmetry improved 
operators as described in appendix [Dj This is to be ex- 
pected because at the level of local operators, increased 
symmetry reduces complications caused by mixing. 



IV. NUMERICAL RESULTS 



A. Mapping out the (I 2 , £P)-plane 




(b) 




FIG. 8: The unrenormalized amplitude A£ nren (7 2 , l-P) ob- 
tained directly from the ratio -R[0^* u _ d [C ; lat ]](P) using the 
sequential propagators with P — (—1,0,0) x 2n/L on the 
coarse- 10 ensemble and applying HYP smearing to the gauge 
fields, (a) real part, (b) imaginary part. 



Following the methods outlined in III B| we have com- 
puted the invariant amplitudes Af nTen (l 2 , l-P) for the 
coarse-10 ensemble, using unrenormalized operators. Ac- 
cording to Table [VT| in the appendix, a straight link cal- 
culation with the operator Or[C/] for T = 74 gives us 
access to A 2 (l 2 ,l-P). Results for A% nTen (l 2 , l-P), in the 
domain we can reach with the available lattice nucleon 
momenta P, are displayed in Fig. [HJ The accessible 
"kinematical" domain is characterized by a triangle with 



an opening angle given by the largest nucleon momentum 
\P\ available in the calculation, see Eq. (24). At I 2 = 0, 



all amplitudes can only be extracted for the single data 
point l-P = 0. The l-P dependence can thus only be 
studied at non- vanishing values of I 2 . Therefore, the x- 
dependence of PDFs cannot be obtained from a direct 
evaluation of Eq. ( fl8| ) on the lattice, in accordance with 
the common knowledge that the lightlike gauge links in 
the gauge invariant definition of PDFs cannot be real- 
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ized on an Euclidean space-time lattice. Nevertheless, 
we will be able to discuss the fcj_-dependence of the low- 
est a;-moment of TMDs, and, beyond that, to draw some 
conclusions about the x-dependence from data at non- 
zero I 2 . 

Coming back to the amplitude in Fig. [HJ we note that 
the real part Re A™ 1011 is dominated by a Gaussian-like 
drop with |/|, while the dependence on l-P at constant \l\ 
features only a slight curvatur e. O ur results for the imag- 
inary part Imi5° ron in Fig. 8b form a surface twisted 
around the |Z|-axis at l-P=0, where the amplitude must 



(a) 



vanish, cf. Eq. (13). The slope of the surface flattens out 
towards larger ]7[7 We will investigate this behavior in 
Section EH 



B. A study of rotational symmetry 

We now study the amplitude A2(l 2 ,l-P) in Fig. 
greater detail for 1° — 0, P — 0. In this case, l-P = 0, i.e., 
the amplitude only depends on the (Euclidean) length of 
the gauge link |Z|. Carrying out the calculation with an 
unrenormalized lattice operator Oy t [C ; lat ], we obtain an 

unrenormalized amplitude A™ 1 ™. Renormalization will 
eventually be based on Eq. (31). However, it is not a 



priori clear to what extent 5m should be independent of 
the direction of the vector I of the link path on the lat- 
tice, since the discretization prescription for the gauge 
link is not (and cannot be) rotationally invariant. Con- 
sider the set of plateau values R[O^ u „ d [C\ at }}(P=0) ob- 
tained from our selection of link paths Cj at . The lattice 
action is invariant under reflections and permutations of 
the lattice axes, i.e., under symmetry transformations 
of the H(4) group. We have checked that the plateau 
values are indeed numerically equal within statistics for 
link paths C\ &t that are equivalent up to reflections and 
permutations of the (spatial) axes. Next, we ask how 
severely continuous rotational symmetry is broken. In 
Fig. [9] we plot the plateau values as a function of the 
quark separation |Z|. To avoid a cluttered plot, we have 
taken the averages over link paths equivalent under H(4) 
transformations. In Fig. [9a] the operator has been evalu- 
ated on the HYP smeared gauge configurations. Here 
the results from step-like link paths and results from 
gauge links on the axes agree very well, and may be de- 
scribed by a smooth, ^-dependent function. A distance 
\l\ where we have a results both from paths along the 
axes and from step-like paths can be found, for exam- 

\l\ = 



pie, at \l\ = 0^42+32 = aV5 2 = 5a = 0.58 fm. We 

find a relative difference of 4 ± 1 percent between the two 
results. 

In the unsmeared case, Fig. |9b| data points from step- 
like links are visibly and systematically lower than data 
points from links along the axes. (At \l\ = 5a, the dis- 
crepancy amounts to 17 ± 2 percent.) We found a very 
similar picture when we studied the breaking of rota- 
tional invariance of the vacuum expectation value of the 



O.N 



0.4 



0.2 



0.0 



X 



X 



X 



X 



X 



x ^x-x-x 



■x.xxxxx- 



0.0 



0.5 



1.0 1.5 

HI (fm) 
(b) 



2.0 



0.S 



O.d 



0.4 



0.2 



0.0 



X 



X 



X 



'X 



x *-x-x-x- 



X-X-X-X.XXX XX- 



0.0 



0.5 



1.0 



1.5 



2.0 



HI (fm) 



FIG. 9: Unrenormalized data obtained for the amplitude 
2A 2yU -d(l 2 ,l-P=0) using the lattice operator 0^[C, lat ] and 
a nucleon momentum P = on the coarse- 10 lattice. Link 
paths coinciding with the lattice axes are marked with a blue 
cross, the red error bars belong to link paths at oblique an- 
gles. The gauge path was constructed (a) on HYP smeared 



gauge configurations, [(b)] on unsmeared gauge configurations. 



gauge link (tr c W lat [C| at ])> on a Landau gauge fixed en- 
semble. As a side remark, we note that a simple cor- 
rection model, the "taxi driver correction" , reduces the 
deviations particularly well in the unsmeared case |36| . 
As a whole, we conclude that rotational symmetry is only 
weakly broken, especially if the gauge link is smeared. We 
rate this as an important indication that the discretized 
operator does indeed approximate the continuum opera- 
tor. In the following, we will analyze nucleon structure 
with the smeared gauge link, and acknowledge a system- 
atic discretization error of the order of four percent as- 
sociated with the violation of rotational symmetry. Last 
but not least, we notice an overall faster drop-off of the 
data with |Z| in the unsmeared case, Fig. 9b than in the 
smeared case Fig. [9a] This can be explained by the fact 
that two different values 8m are needed to renormalize 
the smeared and unsmeared case. 
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Link renormalization 



Method 



In lattice QCD, we work in a cutoff scheme that de- 
pends on the lattice action, with a UV cutoff of the order 
of 1/a. In order to be able to present results for am- 
plitudes that have a well-defined continuum limit, and 
that are independent of the lattice spacing and action, 
we need to renormalize our operator, in particular with 
respect to the self-energy of the gauge link, as discussed 
in Section |III D| The crucial question is how to deter- 
mine 5m in Equations (30 1 and (31). Since we observe 



approximate rotational invariancc for our operator on the 
smeared lattices, we can restrict ourselves to the determi- 
nation of 5m for straight gauge links along one of the lat- 
tice axes. The renormalization of the Wilson line on the 
lattice has a long history in the context of heavy quark 
propagators, where it has been found that the respec- 
tive power divergence requires a non-perturbative sub- 
traction |61j . Calculations in lattice perturbation theory 
[62H64] confirm that the gauge link can be renormalized 
by a factor exp(— 5m L), but will not serve us here to 
determine an accurate value for 5m. Instead, we turn 
to non-perturbative methods. We choose a gauge invari- 
ant procedure based on the static quark potential that 
has been applied in the literature for the renormalization 
of the Polyakov loop [5511681 . Here we outline the basic 
idea. Implementation details are given in appendix [F] 
The static potential V(R) for a system of a heavy quark 
and antiquark with relative distance R can be obtained 
from the asymptotic behavior of the expectation value of 
a rectangular Wilson loop W(R,T) 



W(R,T) = c(R)e 



-V{R)T 



higher excitations , (32) 



where the contributions from higher excitations are ex- 
ponentially suppressed for large T. The Wilson loop is 
renormalized according to 



W Ten (R,T) 



_ e -<5m (2_R+2T)-4;y(90 o ) 



W(R,T), (33) 



where v(90°) is the renormalization constant correspond- 
ing to the 90° corners of the loop. Inserting this form 
into Eq. ( 32 1 shows that the renormalized static quark 
potential 



V rcn (R) = V(R) + 2 5m 



(34) 



obtains a constant offset caused by the self-energy of the 
gauge links in T-direction. Note that we must ensure that 
the loop's gauge links in T-direction are implemented the 
same way as those we use as part of our non-local oper- 
ator. Smearing of the gauge configurations, for example, 
affects 5m. 

A simple renormalization condition that fixes 5m 
would be to demand V ren (i?o) = at some i?n, which 
has to have a fixed value in physical units, see, e.g., Ref. 
[551 155] . An alternative idea [57J EH] makes use of the 



rhu,d/rh 3 


ensemble 


Tmin 


—5m 


i n 


coarse- 10 


smeared 


« 
U 




0.6 


coarse-06 


smeared 


6 


0.1491(31) 


0.4 


extracoarse-04 smeared 


4 


0.1043(94) 


0.4 


coarse-04 


smeared 


6 


0.1554(45) 


0.4 


fine-04 


smeared 


8 


0.1639(35) 


0.4 


superfine-04 


smeared 


10 


0.1578(17) 


1.0 


coarse- 10 


5 


4239('89 1 l 


0.4 


extracoarse-04 




3 


0.361(60) 


0.4 


coarse-04 




4 


0.397(35) 


0.4 


fine-04 




4 


0.382(10) 


0.4 


superfine-04 




5 


0.361(11) 



TABLE II: Renormalization constant 5m from the static 
quark potential. Errors in brackets are statistical. 



fact that the lattice data is quite well approximated by 
the string potential [55] 



string 



(R) = aR- 



TV 

12R 



(35) 



for not too small quark distances R. Matching this form 
to lattice data 8 and demanding C rcn = fixes 5m and 
avoids introduction of another dimensionful constant. By 
setting (7 rcn = 0, we have introduced a renormalization 
condition. In simple terms, it can be understood as the 
asymptotic condition V Tcn (R) — oR — > for large R. 

Applying renormalization with 5m obtained in this 
way, we eliminate the lattice cutoff dependence of our 
gauge links in favor of a reproducible, non-perturbative 
renormalization condition. A future challenge is to find 
the connection of our renormalization condition with the 
scale dependence of TMDs, see also the discussion at the 
end of section IV Bl 



2. Numerical results 

Table |n] lists our numerical results for 5m = 5m/a 
based on matching to the string potential with C rcn = 0. 
We have fit the exponential form Eq. ( 32 1 to Wilson 



loops, where the minimal temporal extent that was taken 
into account is given by T m j n (in lattice units). Most im- 
portant for the following analysis of the invariant ampli- 
tudes are the smeared coarse lattices, where the full set 
of available gauge configurations has been used. The cor- 
responding numbers are shown in bold letters. The other 
lattices serve us to convince ourselves that the method 
works, but do not enter our results on TMDs and could be 
improved with full statistics and larger values of T m j n . In 
particular, the extracoarse-04 lattice may exhibit strong 



8 introducing only a weak dependence on a matching point, chosen 
here to be 1.5ro, in terms of the Sommer scale tq |70| 
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FIG. 10: Renormalized potential for the four smeared lattices 
with m u ,d — 0.4m s . 



discretization errors, and the rather low values 5m ob- 
tained with T m i n = 3,4 may not be reliable. Note that 
our values 8m correspond to — C(/3)/2 in the notation of 
Ref. ESI. 



Figure 10 displays the renormalized potential for four 
lattices with different lattice spacings a but equal ratio 
of quark masses m U( i/m s = 0.4. The data points have 
been corrected for known discretization errors by adding 
\{V^ t (r) — 1/R) as in Eq. ( |F3"| ) in the appendix, and 



the solid lines have been obtained from the model func- 
tion V(R) in that same equation. The curved dashed 
line shows the string potential Eq. (35), plotted for an 
average a. The vertical dashed line indicates the match- 
ing point. The string potential approaches asymptoti- 
cally a straight line through the origin, which we show 
as a straight dashed line in the figure. We see that the 
method yields a renormalized potential that agrees on 
several ensembles of very different lattice spacings. 
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FIG. 11: (a) Yune(-R) evaluated on the smeared gauge config- 
urations for the four lattices with m u ,d = 0.4m s . 
(b) Y l l™(R), renormalized using 8m determined from the 
static quark potential. The gray background highlights a re- 
gion of link lengths R in which lattice cutoff effects lead to 
visible discrepancies between the different ensembles. 



3. Cross-check with open gauge links 



To convince ourselves that the renormalization con- 
stant 8m obtained from the static quark potential renor- 
malizes straight gauge links in general, we study expecta- 
tion values of straight gauge links on Landau gauge-fixed 
ensembles. A convenient quantity to analyze is 

1 (tr c W lat [C r ]) T , 
Y linc (R) = - -In . w iatr C ,i\ ' ( 36 ) 

\ c L 'J/ Landau-gaugc 

where Cy and Ci are straight link paths of lengths R+a/2 
and R — a/2, respectively. Note that the expectation val- 
ues of open gauge links are not meaningful quantities 
without gauge fixing. The renormalization constants Z z 
cancel in the ratio of gauge links, so that the renormal- 
ized quantity is ^[^(-R) = ^Une(-R) + 8m. Indeed, un- 
renormalized lattice results for Yu ne (R) at different lat- 
tice spacings exhibit visible offsets, see Fig. |lla| It is 
encouraging to see that the offsets nearly disappear in 



Fig. |llb| where we have renormalized with the values 5m 
determined from the static quark potential. Except in a 
region roughly below R < 0.25 fm, we find in fact a very 
reasonable agreement of the lattice results for Y^^{R) 
between the different ensembles. We conclude that lat- 
tice cutoff effects become strong for gauge links shorter 
than about three lattice spacings. Therefore, in the fol- 
lowing, we will exclude data points with R < 0.25 fm 
from our analysis. 

A quantitative comparison of Y\i nc (R) at different 
lengths R, different lattice spacings a and an ex- 
trapolation to the continuum a = can provide a rough 
estimate of the size of discretization errors. We perform 
such an extrapolation in appendix[G] The resulting num- 
ber A[(5m]dis = 0.0194 for the coarse-04 ensemble can be 
effectively treated as an uncertainty in the renormaliza- 
tion constant 8m. 
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THE LOWEST i-MOMENT OF TMDS WITH 
STRAIGHT GAUGE LINKS 

A. The ^-integrated correlator and TMDs 



We already stated in section IV A that the restriction 
to the triangle shaped domain in the (|Z|, i-P)-plane given 
in Eq. ( 24 ) precludes us from performing the full Fourier 
transform Eq. (15). However, within our approach, we 



do have access to the x- integral of the correlator Eq. 
i.e., to the lowest x-moment : 

f dx ^\x,k x ;P,S) = J ^e ll± - k± 



(14), 



x ^m(i,p,s) 



i+=r 



=o ' 

(37) 



The above correlator can be parametrized in terms of the 
lowest x-moments of TMDs, cf. Eqns. ^ to (11 ). As an 
example, consider the case of /i, where we define (see 
also Eq. (pj) 



dx fi{x,ki) =2M A 2 (38) 



where /i is the anti-quark TMD defined with respect to 

$ c . Analogously, g^, h±, and hf^ are differences of 

quark- and anti-quark TMDs. On the other hand, fij^\ 

9i> ^Hl^ an< ^ are the sum of quark and anti-quark 
TMDs. 



B. Gaussian fits and renormalized data 



To be able to perform the Fourier transforms Eq. ( 39 ) 
and to renormalize our amplitudes according to Eq. (31 ), 



we follow a simple scheme (This approach circumvents 
potential problems with divergences of the amplitudes at 
|/| = in the continuum limit, see section VD Limita- 
tions of our approach will be discussed later)! 

1. We multiply our unrenormalized data Af nlen (l 2 , 0) 
by the length dependent renormalization factor 
exjp(—5m\l\), using the renormalization constant 
from Table P 9 . 

2. We parametrize the resulting data points in terms 
of Gaussian functions, 



with 



x ^unren 



(Z 2 ,0) 



ill 



_ unrcn 
2 C i,q 



(42) 



A; 




|)^(/ 2 ,0). (39) 



Expressions for the lowest x-moments of other TMDs are 
obtained analogously, in accordance with Eq. (16). Lat- 
tice data for the amplitudes at l-P = are available, e.g., 
from simulations with the nucleon at rest on the lattice, 
P = 0. 

The cc-integral in Eq. ( [37] ) is taken over the whole 
support of $l r ^ (x, fcj_; P, S). The contributions from the 
integration region with x < can be related to anti- 
quark distributions using the correlator $ c defined with 
charge conjugated fields, see Ref. [S] and relation Eq. 
(El ) in the appendix. For straight link paths C sW as well 
as staple shaped gauge links C^ v \ we can decompose the 
x-integrated correlator as 

J dx <5> [r] (x,k ± ;P,S;C) = J^dx $ [r] (x, k ± ; P, S;C) 

+ f dx <S> c ^ c \x,~k ± ;P,S;C), 
Jo 

(40) 



where the parameters c™ rcn , a^ q are obtained from 
fits to the lattice data points. In the fit, we only 
include data points with |Z| > 0.25 fm, to avoid sen- 
sitivity to lattice cutoff effects. It turns out that 
the Gaussian ansatz fits our data reasonably well 
in this range. An exception is the amplitude Ai, 
which appears at subleading twist only. 

3. We determine the multiplicative renormalization 
constant Z^ z by demanding that 



dx / d 2 k ± hjx,kl) = 2A 2 , 9 (0,0) =g v , q 



' <! ' 



where n q is the number of valence quarks (quarks 
minus anti-quarks) . After substitution of the renor- 
malized fit expression for 2^(0,0), the equation 
above reads gv — Z^ el} 1 "™ ■ Since the isovector 
channel is free of contributions from disconnected 
diagrams, we fix Z^ z numerically by setting 



^—1 U — 

^,2 ' „unrcn 



(43) 



-7°7 2 r T 7 2 7°. For T = l, 



and 7 5 , one 



A " y 7 5 , the sign 



where r c = 

finds, r c = r, while for r = 7 M and za 
changes, T c = —T. For the lowest x-moment of TMDs, 
this translates into, e.g., 



dxh(x,kl)- / dxh(x,ki_), (41) 



9 We remind the reader that the renormalization procedure in- 
volves a renormalization condition. In our case, we have chosen 
a condition based on the static quark potential. Changing this 
condition would modify the renormalized data for the amplitudes 
significantly. 
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where n u _<j = 1 and where c™! 1 ^ is directly deter- 
mined froina Gaussian fit to data for the isovector 
amplitude A^u-d- 



4. The renormalization constant Z^,\ 



thus extracted 



from the long-range behavior of A 2 , u -d is applied to 
all amplitudes: We obtain renormalized data points 
from 



Ai, q (l 2 ,0) 



7 -l -Sn 



|/| junrcn 



(Z 2 ,0), (44) 



as well as renormalized fit functions 
1 



i"Gauss 



(I'D 



(45) 



with 



_ 7— 1 „unrcn 
-i,q — ^iS,zH,q 



The prescription above is designed to provide lattice 
scheme and lattice spacing independent results for the 
long-range behavior of the amplitudes Ai. Qualitatively, 
the large- |Z|-behavior of our amplitudes is linked by a 
Fourier-transform to the sm all-| fcj_ |-behavior of the cor- 
responding TMDs, cf. Eq. (H5l to (17 1, since we can 



successfully fit (most of) our data with Gaussians for 
> 0.25 fm, we expect to obtain a reasonable descrip- 



fei 



< 



tion of the corresponding TMDs at small |fcj_ 
1/0.25 fm« 0.8 GeV. 

By restricting the fit to \l\ > 0.25 fm and using 
(smooth) Gaussians to bridge the gap between \l\ = 
0.25 fm and |/| = 0, we effectively regularize any potential 
continuum divergence at |Z| = 0, albeit in a parametriza- 
tion dependent way. This will be important for the def- 
inition and interpretation of (fcj_) "-weighted integrals of 
the TMDs below in section ITCl 

We now discuss results for the coarse-04 ensemble, with 
a pion mass of about 500 MeV. In Figs. 12 and 13 the 
open data points show the unrenormalized amplitudes 
obtained at l-P — 0. From the Gaussian fit to A 2 , we 
determine = 0.938 ± 0.005 sta t. ± 0.042 A[(5A] , where 
the second error is associated with the combined uncer- 
tainly A[<5m] that will be specified in the paragraph be- 
low. The fully renormalized data points are shown as 
solid symbols in Figs. 



12 and 13 The curves and error 
bands correspond to the Gaussian fits after renormaliza- 
tion with Z^ 1 . Data points inside the gray shaded area 
below 0.25 fm have been excluded from the fits. The 
uncertainty obtained from A[<5m] in Eq. (46 1 (see the 



following paragraph) is given by the shaded horizontal 
bands. The fit parameters obtained for the various am- 
plitudes are listed in Table [TlT] Most importantly, we find 
clearly non-zero signals for all amplitudes, even at larger 
distances, except for As and An. Furthermore, the lat- 
tice data points show a high degree of consistency within 
the (in many cases encouragingly small) statistical and 
systematic uncertainties. These results already point to- 
wards rather non-trivial correlations between momentum 
and spin degrees of freedom inside the nucleon. In case 
of the "unpolarized" amplitude A2.U, our data have very 



M 


a 


Oi (fm) 


A 2 , u 


2.0186 ±0.0063 ±0.0008 


1.001 ±0.010 ±0.068 


A 2 , d 


1.0171 ±0.0064 ±0.0005 


0.975 ±0.012 ±0.063 


A 2 ,u-d 


1.0000 


1.029 ±0.018 ±0.073 


A 3 , u 


-0.0978 ±0.0047 ±0.0024 


1.136 ±0.032 ±0.066 


A 3 , d 


-0.0375 ±0.0026 ±0.0009 


1.159 ±0.047 ±0.071 


Az, u -d 


-0.0599 ±0.0037 ±0.0014 


1.125 ±0.044 ±0.065 


A& : u 


-0.9080 ±0.035 ±0.015 


1.207 ±0.036 ±0.089 


A~6,d 


0.2870 ±0.019 ±0.0033 


1.023 ±0.048 ±0.059 


Ae :U —d 


-1.1920± 0.037 ±0.019 


1.164 ±0.026 ±0.080 


A~7,u 


-0.1041 ±0.0064 ±0.0021 


1.151 ±0.047 ±0.074 


A 7 , d 


0.0232 ±0.0038 ±0.0004 


1.079 ±0.12 ±0.063 


A-j tU -d 


-0.1278 ±0.0063 ±0.0025 


1.140 ±0.037 ±0.073 


Ag, u 


-0.0164 ± 0.0048 ± 0.0001 


0.359 ±0.058 ±0.004 


As tU -d 


-0.0178 ±0.0035 ±0.0001 


0.433 ±0.047 ±0.007 


Aq, u 


-0.9268 ±0.030 ±0.011 


1.101 ±0.028 ±0.073 


A~9,d 


0.2636 ±0.016 ±0.0027 


1.057 ± 0.051 ±0.066 


Ag iU -d 


-1.1944 ±0.034 ±0.015 


1.089 ±0.023 ±0.070 


A\Q. U 


0.0881 ±0.0052 ±0.0020 


1.134 ±0.036 ±0.067 


Al0,d 


-0.0137 ± 0.0031 ±0.0003 


1.188 ±0.18 ±0.076 


Aio r u -d 


0.1024 ±0.0054 ±0.0024 


1.139 ±0.033 ±0.067 


Au t u 


-0.0047± 0.0016 ±0.0002 


0.986 ±0.16 ±0.041 


All, u-d 


-0.0045 ±0.0015 ±0.0002 


1.102 ±0.19 ±0.053 




-0.9110±0.032 ±0.0053 


1.058 ±0.035 ±0.072 


^9m,d 


0.2683±0.017 ±0.0015 


1.013 ±0.062 ±0.064 


Ag miU — d 


-1.1822 ±0.034 ±0.0077 


1.046 ±0.027 ±0.069 


A 2 +6,u 


1.1206 ±0.035 ±0.0054 


0.851 ±0.021 ±0.039 


A 2 +6,d 


1.2962 ±0.021 ±0.0088 


0.989 ±0.015 ±0.058 


A 2 +e,u-d 


-0.2451 ±0.034 ±0.0064 


1.622±0.18 ±0.17 


A 2 -e,u 


2.8989 ±0.035 ±0.023 


1.066 ±0.014 ±0.071 


A 2 -e,d 


0.7265 ±0.020 ±0.0041 


0.956 ±0.025 ±0.054 


A 2 -6,u-d 


2.1756 ±0.036 ±0.022 


1.104 ±0.019 ±0.075 




1.0969 ±0.032 ±0.0031 


0.956 ±0.029 ±0.058 


^2+9m,d 


1.2805 ±0.019 ±0.0039 


0.986 ±0.017 ±0.062 


^42+9m,u-d 


-0.1980 ±0.034 ±0.0020 


1.068 ±0.13 ±0.066 


A 2 -gm tU 


2.9113±0.032 ±0.011 


1.024 ±0.015 ±0.069 


A 2 -Qm,d 


0.7483 ±0.018 ±0.0015 


0.958 ±0.029 ±0.059 


^2-9m,u-d 


2.1673 ±0.034 ±0.011 


1.044 ±0.019 ±0.071 



TABLE III: Results from Gaussian fits on the coarse-04 en- 
semble at m v ~ 500 MeV. The first error is statistical. The 
second error includes the statistical uncertainty in Sm and an 
estimate of discretization uncertainties, as given in Eq. ( 46 1 . 



The values for u — d-quarks have been obtained directly from 
Gaussian fits to the u — d data. Note that we have performed 
the conversion to physical units using the values for the lattice 
spacing a given in Table [I] See also footnotes 6 and 7. 



small statistical errors, and we obtain a comparatively 



16 



0.0 0.5 1.0 1.5 0.5 1.0 1.5 2.0 




0.0 0.5 1.0 1.5 0.5 1.0 1.5 2.0 

|1| (fm) |1| (fm) 



FIG. 12: Amplitudes on the coarse-04 ensemble at ~ 500 MeV. We show the unrenormalized data (open symbols), 
renormalized data (full symbols) and Gaussian fits. The uncertainties combined in A[<5m], are given by the shaded horizontal 
bands. 
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FIG. 13: Amplitudes on the coarse-04 ensemble, continued. For convenience, we have introduced a combined amplitude A% m , 
which is associated with the TMD hi. The uncertainties combined in A [5m], are given by the shaded horizontal bands. 
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large value of 3.9 for \ 2 P er degree of freedom. 10 In 
a fit that excludes step-like link paths, x 2 /P-d.o.f is re- 
duced to 2.0, indicating that the the small violation of 
rotational symmetry present in our calculation is to a 
large degree responsible for the high x 2 -value. In the case 
of the twist-4 amplitude Ai Ul we obtain an even larger 
value, x 2 /P-d.o.f = 4.8 both with and without step-like 
paths. In contrast to the case of A% u , the data points 
visually follow a different curve that deviates from the 
Gaussian fit function. The same is true for A\^- We 
conclude that the Gaussian model does not adequately 
describe amplitude A\. Statistical fluctuations are still 
too large to obtain stable fits to Ag^ and An |( j. The 
meaning of the amplitudes A 2 ±e and A 2 ±g m will be dis- 
cussed in section IV E[ 

In order to get an estimate for systematic errors, we 
combine the statistical error A[(5m] sta t and the estimate 
of discretization uncertainties A[<5m]di S of appendix [G[ 



A\Sf 



A[£m] s 2 tat 



A [5m 1 2 



dis 



(46) 



and find that A[5m]j is dominates. It turns out that 
A [5m] mainly affects the widths of the renormalized 
Gaussians, not so much the renormalized c^ q , because 
largely cancel in the process of 



variations in the c"" rcn 



renormalization with . Next, we estimate discretiza- 
tion errors associated with the breaking of rotational in- 
variance. We compare two different Gaussian fits to the 
self-energy-renormalized data for A 2lU - In one fit, we use 
all the data points above |Z| > 0.25 fm, in another fit we 
restrict ourselves to data points from straight link paths 
on the axes. On the coarse-04 ensemble, the relative dif- 
ference in cij 1 ^™ is just 0.6%, and the relative difference 
in cr 2 ,u is 1-5%. Analogous to the case of A[5m], the ef- 
fect on the renormalized parameters c^g is expected to 
be even smaller. We assume that our estimate is also 
valid for the other amplitudes, where it is more diffi- 
cult to make such a comparison due to larger statistical 
errors. In the following, we do not show uncertainties 
from violation of rotational invariance, because they are 
negligible compared to statistical uncertainties and un- 
certainties accounted for in A [5m]. Quantities given in 
physical units are also affected by the uncertainty in the 
lattice spacing a, which is not included in the errors we 
quote. It can, however, easily be obtained by adding 
a relative uncertainty of \d\Aa/a to any quantity given 
in units GeV d or fin . Other sources of errors we do not 
treat here include contributions from excited states in the 
three-point function and the static quark potential, con- 
tributions from disconnected diagrams, and effects of the 
finite lattice volume. Finally, in order to obtain results at 



10 Strictly speaking, we cannot make strong probabilistic arguments 
based on our values of x 2 /p-d.o.f, b ecau se we do not treat po- 
tential correlations explicitly in Eq. ||29[l. 



the physical point, the lattice results as functions of the 
pion mass have to be extrapolated to mjj hys . Although 
we have already performed some preliminary studies with 
respect to the above mentioned issues, they are beyond 
the scope of this initial investigation and will have to be 
left for future work. 

A remaining challenge within our procedure is to asso- 
ciate a renormalization scale with the self-energy renor- 
malization condition we employ. Especially the widths of 
our amplitudes and of the resulting x-integrated TMDs 
are very sensitive to dm, and thus to the employed renor- 
malization condition. We remark that the issue of gauge 
link self-energy appears for any link geometry that con- 
tains space-like sections. Of great interest for future lat- 
tice studies in particular is the development of theoreti- 
cally more accurate definitions of the correlator Eq. ([I]) 
as discussed in the introduction [I] For our purposes, it 
would be important to have subtraction and/or soft fac- 
tors included that cancel the gauge-link self-energies right 
from the start, as discussed already in Ref. |32| . 



C. Interpretation of the lattice results in terms of 
transverse momentum dependent distributions and 
quark densities 



Using Eqns. (38), (39) and analogous Fourier-trans- 



forms for the other TMDs, we can now determine x- 
integrated TMDs from the Gaussian fits to the ampli- 
tudes discussed in the previous section. As an exam- 
ple, for the unpolarized distribution fi we obtain from 



Eqns. (381, (39) 



C9 On 

47T 



(47) 



The result for up-quarks is shown in Fig. |14| Using the x- 
integral of Eq. (16), it is easy to express all x- integrated 
TMDs in terms of the parameters Cj , cr; provided in Ta- 
ble [ill] Note that we have chosen to determine cg m , crg m 
directly from Gaussian fits to the combined amplitude 
Ag m . This way, all resulting expressions for the lead- 
ing twist TMDs are again single Gaussians of the form 
cexp(— k 2 ± /a 2 ). For convenience, we list the numerical 
results for c and a in Table HVl 



In most cases, the widths a turn out to be fairly simi- 
lar. Correspondingly, flavor ratios /j^fcxV/idOkj.) an d 



15a 



and 



15c 



respec- 



orti . In contrast , 



4'i(fcl)M'k fc i) shown in Fes- 
tively, are relatively flat functions'" 

the width of g^\ is significantly lower than that of g^\, 
resulting in a clearly visible slope of the flavor ratio in 
Fig. |15b[ By and large, it is interesting to see that the 
fc_i_-distribution for the down-quarks appear in all three 
cases to be broader than for the up-quarks. In quali- 
tative agreement with our findings, experimental results 
by the CLAS collaboration [71] analyzed using the ap- 
proach of Ref. [75] favor a reduced width of g 1 as com- 
pared to j\. Note that the plots also show results ob- 
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TABLE IV: Numerical results for rr-integrated leading twist TMDs parametrized in terms of Gaussians of the form 
cexp(— fc^/o 2 ), for a pion mass of ~ 500 MeV, straight gauge links, and a renormalization condition based on the static 
quark potential. We also include results for linear combinations of TMDs corresponding to an alternative Gaussian parametriza- 
tion, see section |YE| The first error is statistical. The second error includes the statistical uncertainty in 8m and an estimate 
of discretization uncertainties, as given in Eq. \A&\ . Note that we have performed the conversion to physical units using the 
values for the lattice spacing a given in Table [I] see also footnote 7. 



taincd for the same quantities with an alternative Gaus- 
sian parametrization which will be discussed in section 

EH 

It is natural to think of TMDs as functions that char- 
acterize probability densities of partons in the nucleon. 
Although the probability interpretation is not rigorous, 
see, e.g., Ref. [5], we provide an interpretation of our 
results in this fashion for the sake of an intuitive pic- 
ture. Transverse momentum dependent quark densities 
are introduced as 

P q (x, k±; A, s ± , A, S±) 

= $[(7++A7 + 7 5 - S 3 ^ +J 7 5 )/2]( 2 . i k ± ;P,S) , (48) 

Here the choice of the matrix T = ,(7 + + A7 + 7 5 — 
sHa +: >~f 5 ) ensures projection on the "good" spinor com- 
ponents O [?D and, simultaneously, on the desired light 
cone quark helicity A and transverse quark polarization 
Sj_ [751 HI] • We introduce the following special cases of 
densities: 



Puu,q = ^ p q (x,k ± ;X,0,A,0) = f hq , (49) 

A,A=±1 



PTU,q= 2_v Pq(x,k±',X,Q,0,S±) 



A=±l 



/l,<3 



SjCjiki . 

JVT,q 



m N 



(50) 



PUT,q= \ 51 Pq( X , k ±'i0, S ±, A,0) 



A=±l 



— — hf 
m N 



(51) 



PLL.q = p q (x, fc+; A, 0, A, 0) = - (j 1; q + \A 9l , q ) , (52) 
PTL,q = P q (x, k±; A, 0, 0, S±) 



= 2 (A, 



A 9lT,q 

m N 



Sjtjiki , 

JlT,q 



m N 



)■ 



(53) 



PLT, q = Pq(x, kj_; 0, sj_, A, 0) 



A 



kj_ ■ s±_ , 

n iL,q 



m N 



Kq 



m N 



(54) 
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FIG. 14: /^'(fe!) for up-quarks obtained using the Gaussian 
parametrization at a pion mass m n ~ 500 MeV. The solid 
curve and the statistical error band in blue have been obtainec 



from a Gaussian fit to the amplitude A2, as shown in Fig. 12 
The gray band on the top indicates uncertainties that can 
effectively be expressed as an error in 5m. The gray region 
at large |fcjj indicates that we qualitatively expect strong 
parametrization dependence to set in at |fc^| > 1/0.25 fm ~ 
0.8 GeV. 



PTT, q = p q (x,k ± ;Q,s ± ,0,S±) = - (jx, q 

+ s ± ■ S ± h hq + Sj(2fcjfc ! ~^ S ^ Si fri 



2m 2 N 



lT,q 



— Ki 



(55) 



where the first and the second index of p indicates the 
nucleon and quark polarization, respectively. 

From the x-moments of amplitudes Ai obtained on the 
lattice, we can construct cc-integrated densities p\^ , and 
decompose them in analogy to Eq. ( 40 1 as 



pW(fe x ;A,s x ,A, S±) 



= J dx p q (x,kj_;X,sj_, A, S±) 

= / dx p q (x, k±; A, s±,A, Sj_) 
Jo 

- / dxpq(x, -fex; -A, sj_, A, S±) 
Jo 



(56) 



where the anti-quark density p q is defined as in Eq. ( 48 ) 
but using the correlator of Eq. (El I in the appendix. 



Here the appearance of minus signs in front of p q and 
A accommodates the sign changes in the Dirac matrix T 
after charge conjugation, i.e., T c — — 5(7+ — A7+7 5 — 
sH<j +: >^f 5 ). We conclude that the cc-integrated densities 
Pg 1 ' are differences of quark densities p q and anti-quark 
densities p q of 
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FIG. 15: Flavor-ratios at a pion mass m n w 500 MeV. The 
solid curve and the statistical error band in blue have been 
obtained from the Gaussian fits displayed in Fig. [12] and 



13 The corresponding errors associated with A [5m] are 
shown as a gray band at the bottom. For the dashed curve 
and the band in orange we have used alternative Gaussian 
parametrizations as discussed in section [V E| The respective 
unc ertainties from A [5m] are shown at the top of each plot. 
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(b) 



/i>i)//r>l 
' l i'L(' 6 x)/' l i'<j(^x) fr° m ^9m (solid) and A2±9m (dashed) 



from A2 (solid) and A2±e (dashed) 
from As (solid) and A2±q (dashed) 
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• opposite transverse momentum — fcj_, 

• opposite light cone helicity —A, 

• same transverse polarization s±. 

Strictly speaking, the densities that are integrated over 
x from —1 to +1 are thus not densities themselves and 
can, at least in principle, become negative. 

With the Gaussian x-moments of TMDs from Table 
|IV| as input, we are in a position to draw plots of the 
x-integrated transverse momentum dependent densities 
of quarks in the nucleon. Two particularly interesting 
and statistically well-determined x-integrated densities 
are p^ T and pl^ L . They feature significant dipole defor- 
mations due to correlations in the transverse spins and 
intrinsic transverse momentum, as can be seen from the 



(a) 



terms proportional to gyr and h\ L in Eqns. ( 53 ) and ( 54 ), 
in combination with^our non-zero results for the relevant 
amplitudes A-j and Aio, see Eq. (16). For corresponding 



density plots and their interpretation, we refer to our pre- 
vious publication Ref. [34]. The dipole deformations can 
be characterized by average transverse momentum shifts 
of the quarks, denoted by (k x ) r pL and (fc x )z,r- These are 
defined by ratios of specific moments in x- and fcj_ of the 
densities, as we will discuss in the following section. 

The density interpretation also guides us in our qual- 
itative understanding of the flavor ratio fi u / fid- Ac- 
cording to Eq. pll), we can decompose this ratio as 



(57) 



fi,u( k l) _ Jo dxfi iV (x, k\) - fl dxfi tU (x,k 2 ± ) 



Jo dx fi,d(x, k ± ) - J dxfi td (x,k ± ) 



where, according to Eq. (49), each of the four terms 



on the right hand side has an interpretation as a k±- 
dependent density of unpolarized quarks/antiquarks. In- 
tegrating numerator and denominator individually with 
respect to k± yields the flavor ratio of valence quarks 
n u /nd = 2 in the proton. If f^\(k\) / f^^k 2 ^) were con- 
stant, we would thus expect to find a value of 2. In- 
deed, our result shown in Fig. |15a| is quite close to 2. 
At low |fej_|, the ratio f^u/ffd is slightly higher than 
2, for large |fej_| it drops below 2. According to the 
equation above, the larger ratio at low |fcjj could be at- 
tributed, for example, to an enhancement of the density 
of up-quarks J dx fi >u (k\) a ^ l° w I^J-Ij to a depletion of 
up-antiquarks f 1 dxfi^ u (k 2 L ) at low |fej_|, or to converse 
effects with regard to the down-flavor densities in the de- 
nominator. The flavor ratio for /j 1 ' + shown in Fig. 
|16a| corresponds to the x-integral of pll for A = A, i.e. 
the density of quarks with the same helicity as the nu- 
cleon, minus an antiquark contribution of opposite helic- 
ity, see Eqns. (52) and (56). In this spin-polarized chan- 



nel, we see a strong excess of the x-integrated up quark 
density as compared to the x-integrated down quark den- 
sity. It is well known that up quarks tend to be aligned 
with the proton helicity, while down quarks exhibit the 
opposite behavior. It is therefore not surprising to find 
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FIG. 16: Flavor ratios of ^-integrated densities at a pion 
mass w 500 MeV. The solid curve and the statistical 
error band in blue have been obtained from the Gaussian 
fits displayed in Fig. [12] and |13| The corresponding errors 
associated with A [6 m] are shown as a gray band at the top. 
For the dashed curve and the band in orange we have used 
alternative Gaussian parametrizations as discussed in section 
|VE| The respective errors from A [5m] axe shown at the 
bot tom of each plot. We show up vs. down ratios 



(a) of /j 1 ' + (ft 1 ' from A2, As (solid) and A2±g (dashed), and 



(b) 



°f fi + h[ from A2, Ag m (solid) and A2±9m (dashed) 



a flavor ratio larger than 2 in this channel. However, it 
is interesting to observe that this effect occurs mainly at 
low transverse momentum, as suggested by the notable 
decline of the flavor ratio with Since the Boer- 

Mulders function hi vanishes in the straight link case, 
the combination fi+hi involving the transversity dis- 
tribution corresponds to the density ptt when s± = S± 
and (fc^-SjJ 2 = k 2 j_/2, i.e., on the lines where k± is at 
an angle of 45° with the transverse spin vectors of pro- 
ton and quark. The flavor ratio for this combination is 
displayed in Fig. |16b[ where we observe a similar but 
somewhat less pronounced effect compared to the longi- 



tudinally polarized case in Fig. 16a 
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D. Combined x-fci-moments of TMDs and 
densities 

In the following, we denote the combined x-k±-mo- 
ments of TMDs as 



j[n](m) _ J 



dx x™- 1 I d 2 k ± 



2m 2 N 



(58) 



and analogously for the other TMDs g\, gix, ■ ■ ■ ■ 

As has already been mentioned before, fc^-integrals 
of TMDs taken over the full range of k± are in general 
not well defined due to their asymptotic fcj_-dependence. 
Perturbative calculations show that, e.g., fi(x,k±) ~ 
l/kj_ for large k±, leading to a logarithmically diver- 
gent fex -integral, see e.g., Ref. [38] . Correspondingly, 
in the continuum, the amplitude 2A 2 (l 2 ,l-P) is expected 
to diverge for \l\ — > 0. The required (systematic) regu- 
larization of these potential divergencies will in general 
introduce a dependence on a regularization scheme and 
parameter, e.g. a UV cut-off scale A. Here, we follow a 
simpler, more practical approach and employ the Gaus- 
sian parametrizations of the amplitudes as discussed in 



Section VB which allowed us to perform the necessary 
extrapolation in \l\ to \l\ = 0, and which in turn lead 
to Gaussian (i.e. exponential) fall-offs of the TMDs as 
k± — » oo. With this provisional Gaussian regulariza- 
tion in mind, we can now define a number of ratios of 
fcj^-moments of TMDs and densities that have clear and 
interesting physical interpretations: 
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(59) 
(60) 



giving the well-known axial vector and tensor charges, 
respectively, and 
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(63) 



The first two are the above mentioned transverse mo- 
mentum shifts for longitudinally polarized quarks in a 
transversely polarized nucleon (TV) and vice-versa. For 
later discussions, we have also introduced the transverse 
momentum shift perpendicular to the transverse nucleon 
spin for unpolarizcd quarks (TU), which is given by the 
Sivers function f^p and thus vanishes for straight gauge 
links. We note that the quantities above can be expressed 
in terms of simple ratios of amplitudes, as shown in Eqns. 
p)]) - ([62]) for the case of straight Wilson lines ("sW"). A 



noteworthy advantage of such ratios of amplitudes com- 
pared to individual amplitudes is that they in general 
need no renormalization with respect to the self-energy 
of the gauge link and the multiplicative renormalization 



factor z in Eq. (|31|), i.e. 



Ml 2 ,...) _ A~(l 2 ,...) 



Mi 2 ,...) 



^unren 



(P,...) 



(64) 



due to cancellations of the factors in the numerator 
and denominator. We have to keep in mind, however, 
that we do not evaluate the amplitudes directly at small 
\l\ < 0.25 fm, but rather use the Gaussian parametriza- 
tions to perform an extrapolation to |/| = . There- 
fore, our results can have a residual dependence on 5m, 
and thus on the employed renormalization condition, i.e. 
Qxen _ q Numerically, it turns out that this dependence 
is weak. It is important to note that apart from g^~ d , 
the tensor charge gx as well as the transverse momentum 
shifts are generically scale and scheme dependent quanti- 
ties, due to to the required regularization of the potential 
singularities at very short distances, i.e. the renormal- 
ization properties of the underlying local operators. At 
this point, we are unfortunately not able to relate our 
simple Gaussian regularization to a standard scheme like 
the MS-scheme at a certain scale /Li. This most likely 
requires a detailed theoretical understanding of the be- 
havior of the lattice amplitudes at small |Z|, which may be 
obtained for example using lattice perturbation theory. 
We plan to address this issue in future works. 

The numerical values for the observables given in 
Eqns. (59) to (62) are listed in Table |v| for different fla- 
vor combinations. We note that the value we obtain for 
the isovector axial vector coupling g u ^ d — g^~ d /gv~ d — 
1.192 ± 0.037 ± 0.019 agrees within statistics with the 
value 1.173 ± 0.029 of Ref. [77], obtained using con- 
ventional, local operators on the same ensemble, and is 
also reasonably close to the experimental result g u ^ d = 
1.2694(28) [78 . Our result for the isovector tensor charge 



u,- 

(62) 5t 



g T /g- 



t — d 



1.182 ± 0.034 ± 0.008 turns out 
to be w 10% larger than the value g^~ d ~ 1.06 ± 0.02 
from Ref. [79 obtained for the same ensemble using lo- 
cal operators 11 . This may be related to the fact that 



in the MS-scheme at fi 



4GeV 2 
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observable 


flavor 


value 


gA/gv 




u 


0.450 ±0.018 ±0.008 


Q A i Qv 


[A 2 ± 6 ] 


u 


0.442 ±0.017 ±0.002 


gA/gv 




d 


-0.282 ±0.018 ±0.004 


Q A 1 Qv 


[A 2±6 ] 


d 


-0.282 ±0.018 ±0.001 


gA/gv 




u-d 


1.192 ±0.037 ±0.019 


Oa 1 Qv 


[A 2 ±e] 


u-d 


1.254 ±0.036 ±0.005 


Qt 1 Qv 




u 


0.451 ±0.016 ±0.003 


Qt i Qv 


L4.9-l-Qm 


u 


0.453 ±0.016 ±0.001 


Qt 1 Qv 




d 


-0.264 ±0.017 ±0.002 


Qt 1 Q\f 

a -i / a V 


L4.9-I-Qm 


d 


-0.262 ±0.017 ±0.001 


Qt / Qv 

a 1 I a* 




u-d 


1.182 ±0.034 ±0.008 


Qt 1 Qv 

a - 1 1 a v 


L4.9-|-Qtti 


u-d 


1.201 ±0.034 ±0.002 


{k x ) T L 




u 


69.7 ± 


4.3 ± 


1.4 MeV 


{k x ) T L 




d 


-30.9 ± 


5.1± 


0.6 MeV 


{k x ) T L 




u-d 


172.8 ± 


8.5 ± 


3.3 MeV 


{k x ) L T 




u 


-59.1 ± 


3.5 ± 


1.4 MeV 


(kx) LT 




d 


18.3 ± 


4.1 ± 


0.4 MeV 


{k x ) LT 




u-d 


-138.5 ± 


7.4 ± 


3.2 MeV 



TABLE V: Numerical results for x-fcj_-moments of TMDs ob- 
tained using the Gaussian amplitudes at a pion mass m n ~ 
500 MeV. We also include results corresponding to an alterna- 
tive Gaussian parametrization based on linear combinations of 
amplitudes, as indicated in square brackets, see section [V E| 
The first error is statistical. The second error includes the 
statistical uncertainty in 8m and an estimate of discretization 
uncertainties, as given in Eq. (46 1. The values for u— d-quarks 
have been obtained directly from Gaussian fits to the u — d 
data. Note that we have performed the conversion to physical 
units using the values for the lattice spacing a given in Table 
[I] see also footnote 7. 



the Gaussian para metri zation of the corresponding am- 
plitude Ag m in Fig. 13 in fact overshoots the lattice data 
points at small values of \l\ ~ 0.25fm by « 7 — 10%, 



12 



in contrast to the case of the amplitude Aq in Fig. 
that gives qa- A more sophisticated parametrization of 
the |Z|-dependency of the lattice data for the amplitudes 
could help to resolve this issue. In any case, we interpret 
the outcome of these comparisons as a first non-trivial, 
successful consistency check of our method. 

As we have already discussed in Ref. [33], the aver- 
age transverse momentum shifts, (k x )xL and (k x ) lt (cf. 
Table [V]) turn out to be sizeable and of opposite sign 
for up- and for down-quarks. Moreover, as has been ob- 
served in Ref. [80 , our values are quite similar to the 
results from a light-cone constituent quark model cal- 
culation [5T]. This is remarkable, not only because the 
quark masses employed in the lattice calculation are still 
unphysically large, but also because possible dependen- 
cies on the UV-cutoff scale have neither been investigated 
by us nor in the model calculation. As discussed earlier, 
these dependencies may be weak in particular for quan- 
tities like (k x )xL an d {k x ) lt that can be expressed as 
ratios of amplitudes. It is also interesting to note that 



the gauge link and its geometry do not enter explicitly 
in the calculation of time-reversal even TMDs within the 
aforementioned constituent quark model. 

Finally, we note that as an alternative to the Gaussian 
approach, it is conceivable to regularize the quantities 
defined in Eqns. (59 1- (62 1 by evaluating the ratio at a 
small but nonzero TTT: 



Mo,o) 

1,(0,0) 



reg 



(65) 



For a direct calculation on the lattice, |^ m i n | would have 
to be chosen large enough compared to the lattice spacing 
a to avoid significant discretization errors. 



E. Parametrization dependence using the Gaussian 
prescription 

The simple Gaussian ansatz for the fcj_-dependence of 
TMDs is very successful at parametrizing experimental 
data [82H85] . It also describes the P-dependence of our 
lattice data for the invariant amplitudes at l-P — quite 
well and enables us to perform the Fourier-transform to 
obtain ^-moments of TMDs in a simple way. However, 
this ansatz clearly introduces additional parametrization 
uncertainties. 

In the following case study of parametrization uncer- 
tainties we compare two different ways to use Gaus- 
sians for the parametrization of our data. Consider re- 
integrated densities of longitudinally polarized quarks in 
the longitudinally polarized nucleon 

p ± [ 1 ](fci)^pW(fe ± ;A=±l,A=±l) 



= \{f\ i] (kl)±^(kl) 



A 2 TA ( 



A 



2=p6 



(66) 



In the previous sections, we have discussed individual 
Gaussian fits to A 2 and Aq. This translates into a Gaus- 
sian parametrization of /j 1 ' and with the help of Eq. 
(16 1. Let us label the corresponding results /j 1 ' [A^ auss ] , 
etc. An alternative is^to fit Gaussians to each of the 
combined amplitudes A 2 ±e = A 2 ± Aq. This translates 
directly into a Gaussian parametrization of p^ 1 ', while 
/j 1 ' and now need to be expressed as linear combina- 
tions of Gaussians. Specifically, we obtain 



/i 1] [!grM) 



p+W(*i) 



C2-6C 2 -6 
47T 



g (2/ CT2 _ 6 )^ _|_ 



47T 



(67) 



Note that a single Gaussian function does not change 
sign. Therefore, the alternative parametrization in terms 
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0.6 0.8 
(GeV) 

FIG. 17: The solid curve and error band in blue give the 
relative difference between two different parametrizations of 
/j 1 ' for up-quarks at a pion mass m v » 500 MeV. The gray 
band at the bottom indicates uncertainties that can effectively 
be expressed as an error in 8m. The gray region at large 
| fe_i_ | indicates the scale where we qualitatively expect strong 
parametrization uncertainties to set in. 



of A 2 ±q is in this sense physically better motivated, since 
the quantities ^^(ftj.) have an interpretation as den- 
sities of longitudinally polarized quarks, and should be 
positive [86j . as long as we ignore the (small) contribu- 
tion from anti-quarks, cf. section V C The Gaussian fits 
to data for A 2 +6 an d ^2-6 are of similar quality as those 

we plot for /j 1 ' the rela- 
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for A 2 and A 6 . In Figure 
tive difference between the two parametrizations, namely 

1 - /l 1I [l 2 3auss ]//l 1I [^2±6 S8 ]> as a function of |fe x |. The 
difference between the two parametrizations stays below 
5% for |fej_| < 0.7 GeV , then it rises to an asymptotic 
value of 100% at large |fcj_|. This picture is compatible 
with our qualitative expectations of large parametriza- 
tion dependence beyond |fej_| > 1/0.25 fm w 0.8 GeV. 



Let us now study the ratio 



W(fei)-p" 



(ki 



i(fci) 



(68) 



gj J (fcj) 

as a function of |fe_i_|. In this quantity, both numera- 
tor and denominator become very small at large |fej_|. 
We plot the result in Figure [lSj again comparing the 
two alternative parametrizations. The two results are in 
agreement for \k±\ < 0.6 GeV, at large \k±\ they devi- 
ate strongly. Asymptotically, the curve that corresponds 
to Gaussian g^ 1 and /j 1 ' tends to zero, because has 
a smaller width. The parametrization does not allow a 
sign change of ff^V/j 1 '- O n the other hand, the result 
obtained with Gaussian p + M and exhibits a sign 

change, and tends to —1, because the Gaussian describ- 



ing p 



+ 1 



has a smaller width, so that p 



(68) 



^ ultimately 
It is impor- 



dominates on the right hand side of Eq 
tant to point out that the strong disagreement between 
the two results at large |fej_| is an unavoidable conse- 
quence of the form of the parametrizations, but does not 
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FIG. 18: g l i ] (k 2 ± ) / f[ 1] (k 2 ± ) for up-quarks obtained at a pion 
mass m,r « 500 MeV from two different parametrizations. 
The solid curve, the statistical error band in blue and the 
error associated with A [Sm] shown in gray at the top corre- 
spond to [lef auss ]//l 1] [A§ aass ], while the dashed curve, the 
error band outlined by the dotted curves and the gray error 
band at the bottom correspond to g [ * ] [A§^ BB ]/f[ 1] [AggH- 
The gray region at large |fej_| indicates that we qualitatively 
expect strong parametrization uncertainties beyond [fej_| > 
1/0.25 fm « 0.8 GeV. 



point towards any inconsistencies of the lattice data. In 
this respect, we would like to stress that the same type of 
parametrization uncertainty will at least in principle also 
affect phenomenological TMD parametrizations based on 
experimental data, which are to this date employing 
mostly Gaussian ansaetze for the fcj_-dependence. In 
summary, we see evidence that the relative parametriza- 
tion uncertainty of the Gaussian ansatz becomes very 
large at large |fe_i_|. It appears likely that a better, QCD- 
motivatcd parametrization of the amplitudes at small |/| 
can improve the situation. 



In Figs. 



15 and 16a of the previous section, we have 
always included the result obtained with the alterna- 
tive parametrization based on ^4^g ss . For /j 1 ' ± 
we can introduce an alternative parametrization in anal- 
ogy to Eq. (66) based on Gaussian fits to linear com- 
binations A2±9m = A2 ± Ag m . As before, this ansatz 
seems to be physically better motivated, since the linear 
combinations correspond to (approximately positive def- 



inite) densities as discussed at the end of section VC 
The two types of parametrizations A2 



16b 



A Gauss 



VS. 



In general, we ob- 



^2±9S are compared in Fig. 
serve a rather small difference between them in the range 
< |fej_| ;$ 0.7 GeV. We also include results for the alter- 
native parametrizations in Tables IV and |v| For gA/gv 
and gx / gv j we find differences between the parametriza- 
tions that are in general of the order of the statistical 
errors. For the fits toti-d data, these differences turn 
out to be larger than for the fits to u and d data. 
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VI. TESTING CORRELATIONS IN x AND k ± 

What can we learn from the combined (^-P, /^-depen- 
dence of our amplitudes Ai(l 2 ,l-P) without taking re- 
course to parametrizations and models? A highly inter- 
esting question is if our lattice results for, e.g., A 2 (l 2 , l-P) 
(at least approximately) "factorize" , 



A 2 (l 2 ,l-P) » A 2 (l 2 ,0) A 2 (l-P), 



(69) 



or in contrast show a distinct correlation in l-P and I 2 . 
This is directly related to a corresponding possible fac- 
torization of the x- and fcj_-dependences of the TMDs, 
e.g. 



A(z,fc ± )«A(x) f[ 1] (kl) /M, 



(70) 



where J\f = J d 2 fe±/j i s a normalization factor. 

Model ansaetze based on this assumption are commonly 
employed in phenomenological applications, typically in 
combination with the Gaussian paramctrization of the 
fc^-dependent part, /[^(k 2 ^) j M = exp(— fc5_//x 2 )/7r/i 2 . 
This approach has been used to parametrize experimen- 
tal data of semi-inclusive scattering experiments, see, e.g. 
Refs. 83, 84], and to include effects of intrinsic ("primor- 
dial") parton momentum in Monte Carlo event genera- 
tors, e.g., in PYTHIA and HERWIG++ [S7H89]. Factor- 
ization in x and k\ is a simplifying assumption lacking 
fundamental theoretical justification. Arguments against 
the validity of this assumption have been found in model 
calculations, e.g., in a chiral quark soliton model [10] and 
in a diquark spectator model [91] . see our discussion be- 
low. 



If one of the Equations ( 69 1 or ( 70 1 were to hold 



exactly, it would imply the other one (assuming well- 
behaved functions and integrals). This can be easily 
seen from Eq. (15), which consists of two independent 
Fourier integrals, establishing correspondences I 2 •<->• k\ 
and l-P ■<-» x. The (I 2 , /-P)-factorization thus translates 
into (a:, fc^-factorization of the Fourier-transformed am- 
plitude, and with the help of equation (16), this directly 
implies (x, fc'jj-factorization of f\. 

Analogous arguments connect hypothetical (x, k]_)- 
factorization of other TMDs in Eq. ([16} with (l 2 ,l-P)- 
factorization of corresponding amplitudes A^. 12 

As a first conclusion we note that (x, k]_ )-factorization 
is obviously not in conflict with Lorentz-invariance per se, 
since the parametrization in terms of amplitudes Ai has 
been worked out in a manifestly Lorentz-covariant frame- 
work. We remark that a factorization assumption of the 
momentum-space amplitudes (as defined in, e.g., [5]) of 
the type Ai(k 2 , k-P) — ai(k 2 )cii(k-P) is not equivalent to 



12 For TMDs given in terms of several amplitudes, the latter would 
have to fulfill additional relations among each other. 



the above equations. As a specific example, the on-shell 
approximation Ai(k 2 ,k-P) — S(k 2 )a,i(k-P) discussed in 
Ref. [52] contradicts exact factorization of fx(x,k ± ). 



To study the possibility of a factorization as in Eq. ( 69 ) 



numerically, it is convenient to introduce a normalized 
amplitude 



A™ m {l\l-P) = 



_ Ml 2 , l-P) 



Mi 2 ,o) 



Ml-P) (71) 



and to test whether it is independent of I 2 . We point 
out that the quantity A' lorm (Z 2 , l-P) is renormalization 
scheme and scale independent for finite values of I , since 
both the self energy of the gauge link and the quark field 
renormalization factors of the respective operators can- 
cel in the ratio. As previously, we discard data for very 
small quark separations, |/| < 0.25 fm, to avoid possible 
lattice cutoff effects. In the following, we work with the 
coarse-06 ensemble at m^ « 600 MeV, where we have 
better statistics than on the coarse-04 ensemble due to 
the heavier quark mass and due to a larger number of 
gauge configurations. To reduce discretization errors, we 
use symmetry improved combinations of operators, as 
explained in appendix [D] The effect of this improvement 
turns out to be particularly important for the double ra- 
tios discussed below. Moreover, we make sure that the 
combination of link paths used in the numerator and the 
denominator of Eq. ( |71[ ) are the same up to transforma- 
tions under the hypercubic group H(4). This ensures that 
Sm is exactly the same for numerator and denominator; 
differences in Sm associated with the detailed pattern of 
the link path at the scale of the lattice spacing cancel in 
the ratio. 

\iM-p). 



Figure 
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shows our lattice results for A!, 10 ™ 1 ' 



In each vertical stripe of the plots we show the data at 
constant values of l-P, which is dimensionlcss in natural 
units and can adopt values that arc multiples of 2Tta/L 
with our lattice method. In each stripe, we display the 
results for all available values of |/| = \J — I 2 in the range 
0.25 fm < \l\ < 1.5 fm, with increasing |Z| from left to 
right. For larger \l-P\ only results at larger |/| are avail- 
able, due to the constraint Eq. (24). The data are dis- 



played as filled rectangles representing the statistical er- 
ror bounds, and are drawn with lighter colors for bigger 
errors. 

We find that the data is surprisingly constant within 
the stripes. Taking into account that the errors are corre- 
lated, no statistically significant non-trivial dependence 
on I 2 can be observed in the these plots. Such a depen- 
dence on I 2 would be in conflict with the factorization 



displayed in Eq. (69 1 and Eq. (70) 



Together with the lattice data, we also display results 
from a diquark spectator model [511 ; using formulae 
and parameters given in Ref. [91]. To this end, we cal- 
culate A 2 ornl (l 2 ,l-P) by performing the inverse Fourier 
transform of the analytic model result for fi(x, k 2 ^) nu- 
merically. The choice of the straight gauge link on the 
lattice might be a concern when comparing to models, 
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lattice data 
. model, \l\ = 
model, |/| = 1 fm 
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FIG. 19: Lattice results for the normalized amplitude A^ ""' 1 " 



(I ,l-P), obtained from the coarse-06 ensemble (m, m 625 Me V) 
with HYP-smeared gauge configurations. Each vertical stripe shows results at constant l-P, with values of \l\ ascending from 
left to right. The solid and dashed curves show A™ 1 ™ 1 as a function of I P as obtained from a spectator diquark model [91] for 
several values of |/|. (a) up quarks, real part, (b) up quarks, imaginary part, (c) down quarks, real part, (d) down quarks, 
imaginary part. 
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FIG. 20: Unity minus the double ratio for an exponential 
ansatz for the (fc^)-dependence of the unpolarized TMD for 
up-quarks (employing the parametrization of the GPD H(x, t) 
of Ref. [Ml)- 



however for time reversal even quantities the model cal- 
culations so far do not explicitly include any gauge links. 
Hence, it is difficult to tell at this moment if and how 
this affects the comparison. We remark, however, that 
the lattice calculation has been performed at an unphys- 
ically large pion mass of about 600 MeV, and has not 
been extrapolated to the physical point so far. Neverthe- 
less, we observe a close similarity of the model curves and 
the trend of lattice data. Interestingly, the model results 
for A!} orm as a function of l-P lie relatively close together 
for |Z| = and \l\ = 1 fm. This means that the model, 
when transformed to (7 2 , /-P)-space, also exhibits an ap- 
proximate compatibility with factorization of A2(l 2 , l-P) 
as in Eq. ( 69 1 , at least in the parameter range where 



FIG. 21: Unity minus the double ratio for the diquark spec- 
tator model calculation of f\{x,k\) for up-quarks [91] . 



lattice data is currently available. For larger values of |Z|, 
a possible deviation from the factorization may become 
more visible. 

In order to see more concretely what we can learn in 
principle about the simultaneous dependence of the lat- 
tice amplitudes on (I 2 , l-P) and possible "violations" of 
the approximate factorization, it is advantageous to de- 
fine a double ratio (of, e.g., the real parts of amplitudes) 



R D (l 2 J-P;l 



2 \ — 



Re A? 



Re A? 



(72) 



where Z^ in is the minimal value of I 2 that is available for 
a given l-P and P in our calculation. Clearly, the double 
ratio is strictly equal to unity in the case that the depen- 
dences on I 2 and l-P factorize. We may therefore use its 
variation from unity, 1 — Rjj, as a quantitative measure of 
a potential " violation" of the naive multiplicative factor- 



ization displayed in Eq. (69 1. Furthermore, a cancellation 



of systematic uncertainties and statistical fluctuations is 
even more likely in Ru than in Af ornl . 

To get an idea about what we might expect for the de- 
viation of the double ratio from unity, we show in Figs. 20 



and 211 — i?£>asa function of \l\ for different values of 
l-P, as obtained for two different model-ansaetze for the 
corresponding unpolarized TMD fi(x, kj_) for up-quarks 
in the proton. For a comparison with the lattice results, 
we have, as before, (numerically) Fourier-transformed the 
model-ansaetze to (I 2 , l-P) -space (neglecting sea quark 
contributions by setting f\(x < 0, fcjj = 0), and then 
constructed the double ratio mimicking the restrictions 
in our lattice calculation, i.e. setting P = 2tt / L{n, 0, 0), 
employing typical lattice distance vectors Z, and ensuring 
that \l-P\ < V^P\P\. 

The curves in Fig. [20] are based on an exponential 
ansatz for the (fc^)-dependence and include correlations 
of x and k\ in the form exp(— f{x)k 2 L ). For defmiteness, 
we have chosen the functional form and parameters ob- 
tained in Ref. [94] for the parametrization of the GPD 
H™(x,t), where we have replaced the squared momen- 
tum transfer t by —k]_. This exponential ansatz has the 
right properties in the framework of GPDs, but is un- 
physical in the case of TMDs, and used here just for 
illustrational purposes, i.e., as an example for the type of 
correlations in x and k\ that would be surprising to see 
in our study. As can be seen from Fig. |20| a non-trivial 
signature of the exponential (GPD-like) ansatz in 1 — Ro 
shows up for P = 2vr/L(2,0,0) (for P = 2tt/L(1, 0, 0), 
1 — Rjj is approximately zero in the accessible range of 
variables), where one finds increasingly negative values 
at larger \l-P\ and \l\. 

A distinctly different signature in 1 — Ro is found 
for the TMD f±(x : k 2 j_) for up-quarks from the diquark- 
spectator model calculation of Ref. [91 . In this case, 
comparatively strong deviations from the factorized case, 
i.e., 1 — Rd = 0, are visible already for \P\ = 2ir/L, which 
are, however, positive and hence opposite in sign com- 
pared to the GPD-like ansatz displayed in Fig. [20} Inter- 
estingly, no such clear signature is visible for the corre- 
sponding down-quark distribution. We suppose that this 
is directly related to the fact that the TMD fi(x, k 2 j_) of 
Ref. [91] for up-quarks has a non-monotonic dependence 
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on kj_ at low x, which in turn can be traced back to con- 
tributions of wave functions with non-zero relative orbital 
angular momentum AL Z = ±1. Such contributions are 
absent in this model for fi(x,kj_) for down-quarks. 

Without going into any details, we note that the TMDs 
obtained in the light-cone quark model calculation of 
Ref. [81] also do not factorize, but that at least fi(x, kj_) 
shows a less distinctive signature with respect to 1 — Rjj 
compared to the diquark-spectator model results dis- 
cussed before. In particular, in the model of Ref. [ST] , 
there is no difference between up- and down-quark dis- 
tributions regarding correlations in x and k\. 

Finally, Fig. [22] displays the lattice results for the In- 
dependence of 1 — Rd for eight different values of \l-P\ 
from 7r/10 to 87r/10. Using lattice data points for \l\ > 
0.2 fm, we have constructed Rd for all accessible values 
of I 2 , l-P and the corresponding l 2 ain . Within statistical 
uncertainties, we observe numerically the expected sym- 
metry in IP o -IP, i.e. R D (l 2 ,l-P) = R D (l 2 , -l-P), cf. 
Eq. (13). For the final results, we average over positive 



and negative values to increase the statistics. 

Interestingly, the central values of the lattice results for 
1 — Rd for up-quarks in Fig. 22 show a trend towards siz- 



able, positive values for increasing I at larger l-P, which 
is compatible with the results for the diquark spectator 
TMD model in Figs. [21] However, within the statisti- 
cal uncertainties, the data points are also still mostly 
consistent with zero. Therefore, at present we cannot 
rule out an at least approximate factorization of the 
l-P-, /^dependences of the amplitudes, and the x-, k\- 
dependences of the corresponding TMDs, respectively. 
As a side remark, we note that corresponding lattice re- 
sults for the down-quarks do not show any specific trend 
of the central values at all. It will be highly interesting to 
repeat this study with increased statistics and for larger 
nucleon momenta with, e.g., \P\ = \/2 x 2n/L, 2 x 2n/L, 
and to see if the trend of the central values, pointing to- 
wards a significant correlation in x and k\ as expected 
from certain TMD model calculations, can be firmly es- 
tablished or rejected. 



VII. OUTLOOK 

One of the most exciting challenges for lattice calcula- 
tions of TMDs is to go beyond the direct, straight gauge 
link between the quark fields. This is clearly necessary for 
an understanding of the physics of eikonal phases in pro- 
cesses that involve transverse momentum. The long-term 
goal is to make contact with experimental measurements 
of semi-inclusive deep inelastic scattering and Drell-Yan 
production. Since these experiments are very challeng- 
ing, progress on the lattice in this direction would be 
even more important. What are the principal limitations 
of such calculations? We will need to create a staple-like 
gauge link that resembles the one in Fig. [3aJ i.e., that 
generically runs in a direction v along (or close to) the 
lightcone to infinity and back. First of all, the extent 



of the staple in u-direction, given by the four-vector ryv, 
will always be finite in any practical lattice calculation 
due to the finite lattice volume. By increasing 77 step 
by step, we may hope to find that the data converges 
to a plateau value, which we might interpret as repre- 
senting the limit r\ — > 00. The idea to define the matrix 
element through the limit ry — 00 has already been men- 
tioned in Ref. |32j . Furthermore, on the lattice, we are 
restricted to gauge link structures that have no tempo- 
ral extent, 1° = v° = 0. At a first glance, this might 
seem to imply that lattice calculations with "realistic" 
gauge links are impossible. However, as in the case of 
straight gauge links, we need to establish the connection 
to TMDs using a frame independent parametrization. As 
discussed in Ref. |16j and appendix[C] with an additional 
w-dependence, we now have to deal with 32 independent 
invariant amplitudes, which can depend on the invariants 
I , l-P, rjv-l, {qv) 2 , and r/v-P. The amplitudes defined 
in the limit r\ — > 00 can only depend on ^-independent 
combinations of these invariants. The direction of v rel- 
ative to the nucleon momentum P is essentially 13 given 
by £ = (2v-P) 2 jv 2 , formed from r\v-P and {qv) 2 . For 
finite v-P, the limit of a lightlike staple direction v is 
characterized by |C| — ► 00. On the other hand, insert- 
ing a spatial lattice vector v, we find that £ is bounded 
by < — C < |2P| 2 , where P is the three- momentum of 
the nucleon on the lattice. So although lightlike staple 
links cannot be realized directly on the lattice, the limit 
\(\ — > 00 can still be approached at least in principle by 
choosing larger and larger lattice nucleon momenta. Im- 
portantly, and as already mentioned in the introduction, 
one approach to regularize rapidity divergences in the 
definition of TMDs is to introduce gauge links that are 
slightly off the light cone right from the start, i.e. with 
v 2 0, and hence a finite £. TMDs defined in such a way 
even follow a known evolution equation in the parameter 
C, see, e.g., Refs. [26j EH [95], which allows to evolve to 
arbitrarily large |£|. Based on the above observations, we 
plan to extend our calculations to include staple-shaped 
Wilson lines with varying staple-extents r\, for different 
values of £ employing a larger number of non-zero lattice 
nucleon momenta. To get into contact with the process- 
related TMDs, we will then attempt to extrapolate the 
lattice results to large rj and large \(\, the latter possibly 
with the help of the above mentioned evolution equa- 
tions. This approach should lead to results which may 
be compared in a meaningful manner with corresponding 
results from experimental and phenomenological TMD- 
studies of, e.g., the Sivers effect. To recapitulate, within 
such a formalism, the calculation of TMDs relevant for 
SIDIS or Drell-Yan processes on the lattice could become 
feasible, at least in principle. In practice, one of the fore- 
seeable technical challenges that one has to face in this 
case are diminishing signal to noise ratios for increasing 



The role of the sign of rjv-P is discussed in appendix [(5] 
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FIG. 22: Lattice results for unity minus the double ratio for the real part of the amplitude A2 for up-quarks, for a pion mass 
of ~ 625 MeV. Note that the non-zero nucleon momentum is P — 2tt/L(— 1, 0, 0). 
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nucleon momenta. Another one is the statistical noise 
created by the long gauge link. Furthermore, at present, 
there are also a number of conceptual details concerning 
renormalization of the matrix elements that need to be 
worked out. As pointed out in Ref. [32] . embedding cer- 
tain soft factors in the definition of the correlator could 
cancel the self-energies of the gauge link in an appro- 
priate way. Even without detailed knowledge about soft 
factors, it might be possible to estimate ratios of certain 
fej^-moments such as those in Eq. ( 59 1-( 63 ) , exploiting 



the cancellation of self-energies on the right hand side 
of Eq. ( 65 ) . Especially the transverse momentum shift 



{ky)TU caused by the Sivers function is a promising and 
prominent candidate to investigate with extended gauge 
links on the lattice. 



C,(0) = I, Cz(l) = 0. 

For an arbitrary four- vector w, we introduce light cone 
coordinates w + = (w° + w 3 )/V2, w~ = (w° — w 3 )/V% 
and the transverse projection w± — (0, w , w 2 , 0), which 
can also be represented as a Euclidean two-component 
vector w± = (w±,W2) = (w 1 ,™ 2 ), w±-w± > 0. The 
basis vectors corresponding to the + and — compo- 
nents shall be denoted n and n, respectively, and ful- 
fill n ■ n = 1. The nucleon moving in z-direction has 
momentum P = P + n + {m 2 N /2P + )n and spin S — 
A(P+/m N )n - A(m N /2P+)n + S x , S 2 = -1. We use 
the convention e 0123 = 1 for the totally antisymmetric 
Levi-Civita symbol, and introduce e,-j = e ^ such that 
£12 = 1. 
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Appendix B: Naive Continuum Limit of the Lattice 
Gauge Link 

In this section, we show that the discretized Wilson 
line, given by a product of link variables as shown in 
Eq. ( 20 1 , approaches the continuum Wilson line Eq. ( Al ) 



in the naive continuum limit. 

Consider a lattice path C\ at = (x {n \ x^) that "ap- 
proximates" a continuous, piecewise smooth path Ci of 
fixed length I, By "approximates" we refer to the fol- 
lowing criterium: The path Ci can be subdivided into 
n sections that connect mutually different, path ordered 



(0) on Ci, such that |yW - x 



0(a) 



points y^ n \ . . ., y 
for all i — 0..n. 

Provided the lattice path is not intersecting with it- 
self ( ce;W ^ a;W for all i ^ j ), n must be of order 
i/a for fixed £, since there are 0{t/a) lattice sites a dis- 
tance of 0(a) away from C;. Thus, n grows as a -1 in 
the continuum limit. For reasons of definitcness, we now 
divide the lattice path Cj at into approximately %Jn sec- 
tions, each section connecting approximately the same 
number of consecutive points afW. Consider one of these 
sections, for example, the section running from a point 
x( m ) to x^°K The number of points in this section is 
m + 1 = 0(y/n). For an individual link variable of this 
section, we write 



Appendix A: Conventions and definitions 

Whenever the four- vector I fulfills I 2 < 0, we shall make 
use of the abbreviation |Z| = \J — I 2 . 

In the continuum, a "gauge link" or "Wilson" line is 
given by the path-ordered exponential 



U[d] = V exp(-igj^ de A„(£) 



= V exp \-ig J dX A(C t (\)) ■ C,(A) J . (Al) 

Here the path is specified by a continuous, piecewise 
differentiable function Ci with derivative Ci and with 



U{x {l \x {l - 1] ) = t+igAx {i) ■ A(x (l) ) + 0(a 2 ) 

= l + ig Ax (i) ■ A(x) + 0(a 2 Vn) , (Bl) 



where Ax^ = x^ ^ — x^ — 0(a) with i = l...m, 
and where we used a Taylor-expansion of the gauge field 
A^x) around x = ^ ^ 



A^x) = Afj,(x) + {x- x) y d u A ll (x) + 
= A^(x) + 0(aV^), 



(B2) 



which holds since \x — x\ = 0(a^Jn). For clarity, we have 
kept \fn explicit in our notation, but keep in mind that 
we could formally replace 0(y/n) by 0(a -1 / 2 ). For the 
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product of m = 0(y/n) link variables we then find 



paths are defined as 



U(x {m \x^ m - 1 '>)---U{x il \x w ) 
= t+ig {xW - x {m) ) ■ A(x) + 0{a 2 n) 
= t+ig (2/ (0) - y {m) ) ■ A(y) + 0(a 2 n) . (B3) 

The corresponding section C, of the continuous 
path Ci, running between y( m > and j/ ', reads in ex- 
panded form 14 



U[C<i m ' a) ] = Pexp \ig j^ mo) {A^y) + 0(aVn)}j 
= l + ig (y {0) - y {m) ) ■ A(y) + 0{a 2 n) . (B4) 



Comparing this with Eq. ( B3 ) , we get 



U(xt m \x<> m -V) ■ ■ ■ U(x^,x^) = U[c[ mfi) ] + 0(a 2 n) . 

(B5) 

Analogous relations hold for the other subsections of the 
lattice path and their continuous counterparts. Forming 
the product of these 0(y/n) subsections, we finally obtain 

W latplat] =U [C l ] + 0{a 2 n 3 ' 2 ) U[C{\ , (B6) 

since formally 0(a 2 n 3 ^ 2 ) = 0{a 1 / 2 ) for fixed length £. 



Appendix C: Properties under symmetry 
transformations 



First, consider a general prescription C for the gauge 
paths. Applying Lorentz transformations (P[A]), parity 
transformation (P) time reversal (T), and complex con- 
jugation (f), we obtain the following relations: 



(CI) 

$ [r] (Z,P,S;C) = $ [7 ° r70l (Z,P,-5;C (p) ) , (C2) 
^(l,P,S;C)Y = $hVrV7 1 l(_j i p i s ;C (T)) _ (C3 ) 

$ [F] {l,P,S;C)Y = $ [7 ° rt70l (-?,P,5;C (t) ) . (C4) 



Here the matrices A and A 1 ^ describe Lorentz transfor- 
mations of vectors x^ L — > A^ v x v and spinors ip — > A 1 , 2 i/j. 

For any Minkowski vector w = (w , w) the space inverted 
vector is defined as w = (w , —w) . The transformed link 



'I 



(A) = AC A -i,(A) 



C, (t) (A) = C_ / (l-A) + L 



(C5) 



For straight gauge links U[Ci] = U[l,0], we get C = 
C (L[A]) = C (P) = C (T) = C (t) 5 Le-) the link prescription C 

is invariant. Equation ( CI I then tells us that the correla- 



tor can be decomposed into Lorentz-covariant structures 
weighted by amplitudes Ai(l 2 ,l-P). Equation (C4) es- 



tablishes the relation Eq. (13) between A* (I 2 , l-P) and 
Ai(l 2 ,— l-P). Further relations derived from Eqns. (C2) 
and (|C3| reduce the number of possible non-zero am- 



plitudes, eventually leading to the parametrization Eq. 



As a side remark, we briefly discuss the case of staple 
shaped gauge links in direction v. The paths transform 
according to 

p(ijti)j(L[A]) _£(Aiju)^ p(»7f)l(P) —Q{rfo)^ 

p(j)t))j(T) _Q{-rfo)^ p(i)«)](t) _ (>(riv) _ ^qq^ 

The dependence of the correlator on the direction v leads 
to the appearance of new amplitudes |16j . in total we 
now have 32. Moreover, the amplitudes now depend on 
the Lorentz-invariants I 2 , l-P, rjv-l, (r/v) 2 and rjv-P. The 
amplitudes Ai in the limit rj — > oo can only depend on 
variables that are |?7|-independent combinations of these 
invariants [36, 96]. To obtain a complete set of such vari- 
ables, we divide the invariants by appropriate powers of 
|?7V-P|, a quantity that remains finite in the limiting case 
V = ±n, P + 3> tojv relevant for the discussion of SIDIS 
or the Drell-Yan process. We can thus write the am- 
plitudes as functions A t (l 2 , l-P, v-l/\v-P\, v-P/\v-P\), 
with C _1 = 



: v /\2v-P\ 2 . Inserting Eq. ( |C6[ ) into Eqns. 
(C1)-(C4), we find that the transformations (f) and (P) 
leave v z and v-P invariant, unlike (T), which changes the 
sign of v-P. Therefore, time reversal (T), rather than 
restricting the number of amplitudes, establishes rela- 
tions between amplitudes ..,+1) and Aj(. 1). 
The amplitude with sgn(v-P) = 1 corresponds to SIDIS, 
the amplitude with sgn(u-P) = —1 describes DY. Some 
amplitudes are independent of sgn(w-P), others switch 
sign. Those latter amplitudes lead to "time-reversal 
odd", process-dependent TMDs like the Sivers function 

/it- 



Appendix D: Symmetry improved operators 



14 where, as before, \y — y\ = 0(a^/n) 



Looking at Eqns. (|C5j) , we see that the symmetry 
transformation of the link prescription features a com- 
mon structure consisting of a backward and a forward 
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r (Eucl.) T (Mink.) 
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f [7i 1 72] 
73 

|[7i,73] 
|[72,73] 
-7475 

74 
I [7i, 74] 
I [72, 74] 

7375 
|[73,74] 
-7275 
7i75 
75 



—47 

-i(ia 03 7 5 ) 

-«7 3 
i(ia 02 7 5 ) 



- J ( ? a 01 7 5 ; 



7°7 5 



ia 23 7 5 



-ia 13 7 5 



i7 3 7 5 



ia 12 7 5 



-i7 2 7 5 



i7 1 7 5 



m N 
E(P) 
i 



E(P) 
m 2 N 



A 2 P 1 + ^A 3 l 1 



E(P) 



E(P) 



Ash 



-1A9 - im%Au (Z3) 2 



E(P) ° 
im 2 N An hh 

—im 2 N An hl 3 

irriN A7 I3 

A 2 

o 

Ain h 



E(P) 

1 



E(P) 



A 9 Pi 



E(P) 



E(P) 



im N 7 *m N ~ 2 

A 6 ~ TTFTvT ^8 (h) 



E(P) 



E(P) 
im% 
EjP) 



E(P) 



Ashh 



m » Ashls-^ArlsP, 



E(P) 



E{P) 







TABLE VI: Plateau values of the ratios R[Or n [Ci]](P) for 
straight gauge links Ci in terms of the amplitudes Ai. Here we 
employ the LHPC conventions for T 2pt = T 3pt = (1 +7 4 )(1 + 
*7s73)/2, i.e. the nucleons are spin-projected along the z- 
axis. We choose the nucleon momentum P — (Pi, 0,0), and 
the quark separation is I — (li, 12,13), h = 0. 



When constructing a discretized version of the gauge 
link operator, we can reduce discretization artefacts by 
preserving those symmetry transformation properties of 
the link path that have a correspondence in discrete 
Euclidean space. In this context, it is convenient to 
represent the discrete link path as a sequence of shifts 
of one lattice unit. Let 1, 2, 3, 4 denote vectors of 
length a along the four lattice axes. A lattice link path 
may thus be represented as Cj at = [s^ n \ . . . , sW], with 
s (i) e {_4 ; ^ _i 5 i 5 . . . ; 4}. The sample link path of Fig. 
[i] is given by C\ at = [1,2,1,1,2,1,1,2,1]. On the lat- 
tice, the Lorentz group and parity are replaced by the 
hypercubic group [371 [55] 

H(4) = {(6,7r) I h, b 2 , 63,646 {0,1}, 7r€S 4 } , (D3) 

where S4 is the set of permutations of {1,2,3,4}. The 
action of a given group element h = (b,Tv) of H(4) on a 
link path is given by 



fl W _> s '(i) 



-(-l) b ^(4) 

(-l) bl ^(l) 
(-1)^7^(4) 



M = 



,(0 



,(0 



(D4) 



i.e., H(4) permutes axis labels and inverts the direction of 
lattice axes. This defines D(h). Hermitian conjugation 
(f) of the matrix element reverses the ordering of the 
shifts and negates them: 



D{\) [s 



(n) 



.,S 



(1)1 _ 



»1 



(D5) 



The representation D is deduced from the transformation 
behavior of I — Y17=i 

Dih- 1 ) l = ((-l) bl /. (1) , . . . , (-l) & ^. (4) ) , (D6) 
I){\)l = -l. (D7) 

The operation f is its own inverse and commutes with 
any h £ H(4). Thus we can define a larger group 



transformation that leaves, as a whole, the vector be- 
tween start and end point of the link invariant. In math- 
ematical terms 



C l ^D{g)C 



Dig- 1 ) I 



(Dl) 



Here g is a group element of one of the respective sym- 
metry groups, i.e., Lorentz-transformations, parity trans- 
formation, time reversal or Hermitian conjugation. The 
representation D(g) of that group element acts on the 
link path, while the representation D(g) acts on vectors. 
The representation D can be deduced from D by look- 
ing at the transformation behavior of the vector between 
start and end point of a link path, i.e., 



D(g)l:= [D(g)d] (0) - [D(g) d] (1) 



(D2) 



G= [J {h,\oh} 

heH(4) 



(D8) 



The function C lat we use to determine the link path for 
a given vector I is a Bresenham-likc algorithm that pro- 
duces a step-like path close to the straight continuum 
line. It turns out that, in general, this alorithm is not 
invariant under transformations of the form Eq. (Dll. 



However, it is simple to form a superposition of gauge 
links that has the desired properties: 



—y 



u 



(D9) 



All gauge links in the above superposition run from I to 0. 
Thanks to the properties of the algorithm C lat , the above 
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sum does not contain link paths that have an extent in 
the Euclidean 4-direction. Performing the substitution 
Eq. 1D1|) on the right hand side of the equation above, 



we obtain 



#G 



gee 



lat 

D(g-i)D(g-i)l^ 



# 



gee 



lat 

D((gc,g)- 1 )l 



#G 



sec 



(D10) 



its direction: 



V exp y dX A(d(X)) ■ Ci(\) 
V exp(+ig [ d\ A{d{\)) ■ Ci{\) 



V exp l^+ig J dX A(Ci(l - A) J • C;(l - A) 
P exp / dA AfCi(A)) -di(A) 



(E3) 



because G o g — G. So is indeed invariant under 



transformations Eq. (Dl) for any g € G 



In practice, the sum of Eq. ( D9 ) contains typically only 



a few distinct link paths. We evaluate three-point func- 
tions for all these different paths. In the final analysis, 
we form the superpositions using appropriate weights for 
the individual paths corresponding to their multiplicities 
in the sum. 



Appendix E: Charge conjugated operator 

In the presence of a general link path C, a gauge invari- 
ant definition of the correlator $ c of Ref. [9] is obtained 
by applying charge conjugation G to the whole operator: 



where Ci (A) = C; (1 — A). Using translation invariance, we 
obtain 

$ c W(l,P,S;C) 

=i(P,5|g(0) HVrW) U[Ci\q(t)\P,S) 

=\ (p, s\ q(-i) (- 7 Vr T 7 V) m - 1} q(0) |P, S) 
= $ c[ - 7 Vr T 7 V] ( _ /; Pj s . C (t) } ; (E4 ) 

because C,(A) — I =C_ ( _ i) (l- A) + (-/) =C m l (X). Carry- 
ing out the Fourier transform with respect to I, we arrive 
at Eq. pil. 



Appendix F: Implementation details of link 
renormalization 



-ik-l 



X-(P,S\ Cq^TUid] q(0) C \P,S) 
= $[-7Vr T 7 Y](_ fci p )S;C (t)) . (ei) 



where the conjugated link path C^) is defined in Eq. 
(C5). The straight gauge link and the staple-shaped 



gauge link turn out to be unaffected by the charge con- 
jugation, c sW <t) = c sW , cww = e<». 

For completeness, we show the proof of the third line 
of the above equation. Using C#(i)C = —A^(x), 
Cq{x)C = i7°7 2 g T (a;), Cq(x)C = q T (x)ij°"f 2 , where t 
is acting on Dirac and color indices only, we get 



G q{l)TU[Ci] q(0) C 

= 9 T (i) 17 Vrii[c,]' n°7 2 5 T (o) 
= - g(o) ( i7 °7 2 r l 7°7 2 ) T u\ctf qii) 



(E2) 



In the last line, we have used that fermion fields anti- 
commute. Denoting reverse path-ordering V, we find 
that the Hermitian conjugate of the gauge link reverses 



We calculate rectangular Wilson loops on the lattice 



W ut (r,T) = l(tr c U lat [C r . T }) 



(Fl) 



for closed paths C r t as depicted in Fig. 23 Here r is a 



spatial vector between lattice sites. For the correspond- 
ing spatial sections of the gauge link, we use step-like 
paths as in Section [III A| For large enough T, 



IU lat (r, T) « c(r) exp (~U lat (r) T) 



(F2) 



Taking lattice data at fixed r and a range of values T 
enables us to determine U lat (r) and c(r) from an expo- 
nential fit. To obtain a smooth interpolating curve of the 
static quark potential as a function of R = |r|, and to 
reduce discretization errors, we follow Refs. [45] [99] and 
fit the functional form 

U lat (r) = aR-a/R + C -X (^ at t (f) - 1/p) (F3) 
= V{R) 

to the data obtained for U lat (r). Here the hat " in- 
dicates that the respective dimensionful quantity is ex- 
pressed in lattice units. The potential U p lat t (r) is ob- 
tained from single gluon exchange between the temporal 
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(a) 



FIG. 23: Rectangular Wilson loop in the calculation of the 
static quark potential. 



links in lattice perturbation theory. The corrective term 
[lOOj proportional to A, associated with breaking of ro- 
tational invariance, becomes negligible for R > 3a. For 
the calculation of Vp^. t (r) we use the inverse gluon prop- 
agator of the MILC action [lOlj . and, if the potential 
is calculated on smeared gauge configurations, the ap- 
propriate HYP smearing coefficients hnp(k) from Ref. 
|102j . Once the fit parameters a, a, C and A have been 
determined, we obtain 5m from equating the renormal- 
ized potential V Ien (R) = V(R) + 2 5m with the string 
potential String (R) = &R — n/12R at a matching point 
R = 1.5 r = 1.5 ^(1.65 -a)/cx : 



25m 



-C 



1 

L~5 



1.65 



7T 

12 



(F4) 



Appendix G: Estimating discretization errors from 
the gauge link 
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A comparison of Yr inc (R) defined in Eq. (36 1 for dif- 
ferent lattice spacings allows us to get an idea about the 
size of discretization errors coming from the gauge link. 
While the renormalized quantity Yj[°"(i?) must be inde- 
pendent of the lattice action, smearing and the lattice 
spacing, 5m and Yy mc (R) can assume different numerical 
values for different lattice spacings: 

= W^S «i) + Mai) = Y lins (R; a 2 ) + 5m{a 2 ) 

Thus the right hand side of the difference 

5m(a 2 ) - 5m(ai) = Y Xirie {R) a x ) - Y line (R;a 2 ) (Gl) 

should be R- independent up to lattice artefacts. We es- 
timate the latter by comparing two different link lengths 
Ri and R 2 : 



A(R 1 ,R 2 ,ai,a 2 ) 



'■ {Y lhlc (Ri;ai) 
(Y Vlne (R 2 ; ai) 



,(i?i;a 2 )) 
3(^2; 02)) 



(G2) 



As a technical note, we mention that we employ a spline 
interpolation in order to be able to evaluate Yu nc (R) at 
arbitrary values R. If there were no discretization errors 
at all, A would be zero for any choice of a\, a 2 , R\ and 
R 2 . We remark that Yn ne (R) can naturally provide an 



FIG. 24: (a) Ri -dependence of A(Ri, R2, ai, 02) for fixed ai 
and i?2- The dashed line corresponds to the superfine-04 en- 
semble with d2 ~ 0.06 fm, the solid line to the fine-04 en- 
semble with ci2 ~ 0.09 fm and the dotted line to the extra- 
coarse-04 ensemble with 02 ~ 0.18fm. |(b)| a2-dependence of 
A(i?i, R2, 01, ai,) for fixed ai, Ri and R%. The solid data 
points correspond to the data points extracted at Ri = 
0.25 fm in the figure above. The curves with statistical er- 
ror bands are fits to A assuming discretization errors ~ a v . 
The data point at a ~ 0.18 fm has been excluded from the fit. 
The data points with crosses indicate the extrapolated values 
A[5m] diB at a 2 = 0. Note: An error of A[Sm] Ais = 0.01 GeV 
corresponds to an uncertainty of about 2% in the width 2/02, u 
of the Gaussian we obtain for the i-integrated unpolarized 
distribution /{^(fex)- 



alternative way to fix 5m, e.g., with a (gauge dependent) 
renormalization condition ^[^"(^0) = for some fixed 
length i?o- This has already been suggested long ago in 
Ref. [1031 HD3). Comparing with Eq. (|GT|), we learn that 
A can be understood as a discrepancy in the values 5m 
needed to renormalize Yu nc at two different link lengths 
Ri and R 2 . Our goal here is to estimate discretization 
errors for the coarse-04 lattice, so we need to compare 
Yi inc determined on the coarse-04 lattice (ai 0.12 fm) 
with the other other -04 ensembles (a 2 « 0.06, 0.09, and 
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0.18fm). We choose i?2 = l-5r = 0.70fm, the same 
length scale we use as a matching point in our determi- 
nation of 5m from the static quark potential. Figure |24a| 
shows A(i?i, i?2=l-5ro, ai=0.12 fm, 02) for the different 
available lattice spacings a 2 as a function of R\ . We find 
that the magnitude of A and its slope are largest when 
Ri is small, i.e., when Ri is of the order of a few lattice 
spacings a± or 02- This finding corresponds to the dis- 
crepancies already observed in Fig. lib in the region R < 
0.25 fm and leads to the conclusion that very short gauge 
links suffer from significant discretization errors. We now 
choose Ri = 0.25 fm, i.e., the shortest length of gauge 
links we accept in our TMD analysis. The correspond- 
ing values A(i?i=0.25fm, R2=1.5ro, ai=0.12fm, a%) give 
rise to the data points with statistical error bars at the 



dashed vertical line on the left in Figure 24a The same 
data points are plotted with respect to 02 in Fig. 24b 
Assuming discretization errors of 0(a p ), we have per- 
formed one-parameter fits of the form 



A(iii,.Ra, 01,03) « c(a 2 p - ai p ) 



(G3) 



to the data points in Fig. |24b| At present, we do not 
know the order of convergence p. Appendix [S] shows that 
p > 1/2 in the naive continuum limit. We have tried out 
fits with p — 1/2, p = 1 and p — 2, always excluding the 
data point from the extracoarse-04 lattice from the fit. In 
order to estimate discretization errors for the coarse-04 
lattice, we use the above fits to extrapolate A to a 2 = : 



A[5 



m\ di . 



lim A (ili. 

a 2 — >0 



#2, 0,1, O2) 



(G4) 



where R\ = 0.25 fm, R2 = 1.5ro and ai s» 0.12 fm are 
kept fixed. We can interpret A[(5m]di S as the size of 
a spurious R-dependence of Sm that appears when we 
match Yiinc at finite lattice spacing a\ to Yn nc m the 
continuum over a range of link lengths between R\ and 
i?2- Thus A[5m]di s can be effectively treated as an uncer- 
tainty in Sm. For the three different values of p, we obtain 
from' the fits A[5m} dis = 0.0573(59) sta t GeV, A[<5TO] dis = 
0.0323(34) stat GeV, and A[Sm] dis = 0.0200(21) stat GeV, 
respectively. For our presentation of numerical results 
in section |Vj we select the value obtained from the as- 
sumption of 0(a) convergence: A[(5m]di s = 0.0323 GeV, 
or A[5m]dis = 0.0194 in lattice units. With respect to 
our analysis based on a Gaussian parametrization, the 
main effect of A[<5m]di S is an additional uncertainty in 
the widths &i >q of the amplitudes A ijq (l 2 , 0). 

We remark that our determination of A[5m]di s is based 
on open gauge links 14[C{\ evaluated on a gauge fixed 
ensemble. Discretization effects of the complete gauge 
invariant operator q(l)YU[Ci]q{Q) might be different, es- 
pecially for short gauge links. Our value A[<5m]di S deter- 
mined with open Wilson lines can thus only serve as an 
order of magnitude estimate of potential discretization 
errors. 



Appendix H: Expansion in terms of local lattice 
operators 

The nonlocal lattice operators studied in this work can 
be written as weighted sums of local operators involv- 
ing higher derivatives. It is well known that due to the 
loss of translational and rotational symmetries on the lat- 
tice in particular, the operators with two or more deriva- 
tives will mix with operators of lower mass dimension 
under renormalization. This type of mixing involves in- 
verse powers of the lattice spacing, and hence the respec- 
tive contributions have to be subtracted explictly before 
the continuum limit can be taken, which is in practice a 
difficult task. The question then naturaly arises if and 
how these observations can be reconciled with the known 
renormalization properties of a manifestly non-local op- 



erator as explained and used in section HID Although 
we are not able in the course of this exploratory study 
to provide a definite answer, we will briefly explore this 
question in the following and at least show that our renor- 
malization prescription of the non-local operator on the 
one hand, and operator-mixing within an expansion in 
terms of local operators on the other, are not in any ap- 
parent contradiction to each other. 

To keep the discussion simple, we consider here a non- 
local operator with a straight-link of length I in the di- 
rection of the unit vector e„ 



O r (ee^ = q(0)TU[0Je^}q(£e^ 



(HI) 



Our discrete representation of Or(^e M ) on the lattice is 

[O r ("A)] lat = 9(0) T U(0, {£)■■■ U{{n - n£) q{nfl) , 

(H2) 

where n — £/a. Together with a discretization prescrip- 
tion for the covariant derivative on the lattice, e.g. 

DJ(x) = -{U(x, x + p)f(x + A) - f(x)} , (H3) 
a 

we can write [Or(^A)] lat as a weighted sum of local lattice 
operators: 

[Or(nA)] lat = g(0)r (a^ + 1)"?(Q) 



E 

k=a 



9(0)r££g(0) . (H4) 



= [0£' fc ] lat 



To simplify the discussion of operator mixing, we only 



consider mixing of operators [Op' fc ] lat among themselves 



[O 



\i , k-\ lat 



3=0 

= a- k Z k0 [Op ]™ + --- + Z kk [0£' k Y ea + 



(H5) 

Here the powers of a required to render the mixing coeffi- 
cients Z k j dimensionless can become negative, the "worst 



3G 



case" being the potential mixing with the derivative-free 
operator [Op'°] lat . Inserting the above expression into 
the second line of Eq. ( H4 1 yields 



9=0 lfc=0 ^ ' J 



Cj(n) 



(H6) 



For the discussion of the continuum limit, it is at this 
point important to distinguish two cases: 

1. keeping n = Ija fixed as a — > 0, 

2. keeping £ fixed as a — > 0, i.e. sending n — > oo. 

In the first case, it is easy to see that within the op- 
erator expansion in Eq. |H6[ inverse powers of a due to 
mixing are not an issue, since they no longer show up 
explicitly. Clearly, in the continuum limit, the phys- 
ical extent £ shrinks to zero, and only the oper ator 



H6 



[Op'°] rcn contributes on the right hand side in Eq 
while [Or(^A)] lat f° r fixed n is just the discrete represen- 
tation of a local continuum operator. This local interpre- 
tation of [Or(^A)] lat is, however, not the one relevant for 
this study. 

We now turn to the second case, where the length £ is 
kept fixed. As a — > 0, n — > oo, the number of terms in 
Eq. (H6) increases, and due to the quickly growing bi- 



nomial coefficients, the coefficients Cj eventually receive 
infinitely large contributions. Without detailed knowl- 
edge about the mixing coefficients Z^j , we cannot derive 
the renormalization properties of the non-local lattice op- 
erator from Eq. (H6). It is essential to realize, however, 



that the renormalized form of the non-local operator is 
known, both in the continuum 55-59!, and on the lattice 
from heavy quark effective theory in the static quark limit 
[6D - K31 ITU3"] . Restating Eq. (31 ), the non-local operator 
can be written in terms of the renormalized operators as 



[Or(nfi)] M = Z 9t ,e T - 



[O r (£e, 



(H7) 



Inserting this into Eq. ( H6 ) , we find that 



[Or(%)] ren = £ Z*i 



1 g— n 5fi 



c j (n = £/a)£ : >[0£ J Y cn 



(H8) 



With linearly independent [Op J ] rcn , and assuming a 
marginal (not power-like) a-dependence of the renormal- 
ized operators, one finds that for fixed £ the coefficients 
Cj have to scale in unison with the lattice spacing a in- 



dependent of j, according to 

Cj cx Zip. z e 



Srh £/c 



(H9) 



Such an exponential scaling of the Cj is indeed not an 
implausible scenario and can be driven by the binomial 
coefficients, cf. Eq. |H6| We conclude that a simple di- 
mensional analysis does not reveal any obvious conflict 
between mixing of local operators and the renormaliza- 
tion properties of our non-local operator. By evaluating 
the non-local operator directly, we apparently bypass the 
severe l/a™-mixing problem that complicates the com- 
putation of individual local operators with higher deriva- 
tives on the lattice. 

As a final side remark, we note that the last line of Eq. 



( H4 ) can be simply rewritten as 



[O r (n/!)] lat = f> fe ( £/ °\l/a)- k q(0)TD^q(0) 



k=0 



(H10) 



which has the form of an operator product expansion 
(OPE) |105l I106| in terms of a complete set of local op- 
erators Oi(0) and dimensionless coefficients Ci{£\) (see, 
e.g., chapter 18.3 of Ref. [TO?] ) 



E 



€*- 3 <5,(*A) [O,(0)] r 



(Hll) 



Here, di denotes the canonical mass dimension of op- 
erator Oi, and all renormalized operators in the above 
equation depend implicitly on the renormalization scale 
A. Unlike an OPE in the continuum, the expansion on 
the lattice Eq. (H10) terminates after a finite number 



of operators, but is nevertheless an exact identity among 
lattice operators. For k -C l/a, the binomial coefficient 
is C[ at (£/a) w 1/fc! such that the first terms in the sum 
remind us of a regular Taylor expansion. 



Interestingly, a strategy proposed to overcome issues 
of operator mixing in the calculation of higher moments 
of structure functions [108 involves lattice correlators 
that are quite similar to those employed in the study 
at hand. This strategy introduces a bi-local operator 
^(1)^^(1) #(0)7 1/ g(0) with a fictitious heavy quark field 
^. The connection to our approach can be seen in the 
static quark limit —> oo, where the field ^ can be in- 
tegrated out and 'J'(O) essentially becomes a Wilson 
line in 4-direction. The strategy of Ref. [108] requires 
a continuum extrapolation and interpretation of the bi- 
local operator before local operators are determined from 
the matching to an OPE, thus avoiding complications re- 
lated to the reduced symmetries of the lattice. 
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